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This  thesis  conducts  an  in-depth  study  of  the  computational  issues  associated  with  solving 
a  set  of  coupled  discrete-time  Riccati  equations.  Briefly,  the  organization  of  this  study  is  as  fol¬ 
lows.  First,  the  problem  is  motivated  by  discussing  two  game  situations  which  give  rise  to  cou¬ 
pled  discrete-time  Riccati  equations.  Next,  the  computational  aspects  of  solving  these  coupled 
equations  are  investigated.  Finally,  algorithms  and  software  are  produced  that  iterate  these 
equations  in  a  numerically  robust  and  computationally  efficient  manner.  The  thesis  carries  the 
coupled  Riccati  problem  from  formulation  to  software  implementation  with  several  theoretical 
advances  along  the  way.  However,  the  major  contribution  of  this  work  is  the  Riccati  solution 
method  -  i.e..  the  algorithms  and  software  which  solve  the  problem.  As  the  algorithms  are  for¬ 
mulated.  structured,  and  subsequently  coded,  the  software  engineering  factors  that  influence 
good  software  design  are  addressed.  Furthermore,  the  coupled  Riccati  software  developed  here 
is  integrated  into  a  well-known  Computer-Aided  Design  (CAD)  software  package.  Thus,  the 
informal  computer  user  has  easy  access  to  software  which  solves  both  single  and  coupled  Ric¬ 
cati  equations. 
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INTRODUCTION 


This  thesis  conducts  an  in-depth  study  of  the  computational  issues  associated  with  solving  a 
set  of  coupled  discrete-time  Riccati  equations.  Briefly,  the  organization  of  this  study  is  as  follows. 
First,  the  problem  is  motivated  by  discussing  two  game  situations  which  give  rise  to  coupled 
discrete-time  Riccati  equations.  Next,  the  computational  aspects  of  solving  these  coupled  equations 
are  investigated.  Finally,  algorithms  and  software  are  produced  that  iterate  these  equations  in  a 
numerically  robust  and  computationally  efficient  manner.  The  thesis  carries  the  coupled  Riccati 
problem  from  formulation  to  software  implementation  with  several  theoretical  advances  along  the 
way.  However,  the  major  contribution  of  this  work  is  the  Riccati  solution  method  -  i.e.,  the 
algorithms  and  software  which  solve  the  problem.  As  the  algorithms  are  formulated,  structured, 
and  subsequently  coded,  the  software  engineering  factors  that  influence  good  software  design  are 
addressed.  Furthermore,  the  coupled  Riccati  software  developed  here  is  integrated  into  a  well- 
known  Computer-Aided  Design  (CAD)  software  package.  Thus,  the  informal  computer  user  has 
easy  access  to  software  which  solves  both  single  and  coupled  Riccati  equations. 

1.1  Motivation 

This  section  motivates  the  computational  study  of  coupled  discrete-time  Riccati  equations. 
Historical  background  is  presented  to  give  the  proper  setting.  A  short  discussion  of  the  scope  of  the 
thesis  contribution  follows  to  provide  some  breadth  to  the  findings.  But  the  set  of  computational 
issues  is  the  key  focus  of  this  work  and  the  software  design  and  implementation  are  of 
fundamental  concern.  Hence,  the  Computer-Aided  Control  System  Design  (CACSD)  field  is 
introduced.  Then,  the  L-A-S  CACSD  language  is  reviewed. 
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1.1.1  Background 


The  theory  of  optimal  control  has  reached  a  significant  level  of  maturity  as  evidenced  by  the 


number  of  textbooks  written  on  the  subject  [e.g..  1.2.3].  A  classic  problem  often  discussed  in 


introductory  courses  is  the  Linear  Quadratic  (LQ)  regulator  using  state  feedback.  Here,  one  finds 


that  the  optimal  control  is  obtained  by  solving  a  single  Riccati  equation.  Originally,  the  solution  to 


the  dual  (filtering)  problem  was  given  by  Kalman  [4].  Subsequently,  both  continuous-time  [1.2.3] 


and  discrete-time  [5.6]  versions  of  the  Riccati  equation  have  been  studied  in  detail.  A  key  feature 


of  optimal  control  problems  is  that  there  is  a  single  control  agent  or  Decision  Maker  (DM). 


A  more  interesting  situation  occurs  when  there  are  two  or  more  DMs  controlling  the 


underlying  dynamic  system.  A  straightforward  generalization  of  the  state-feedback  regulator 


problem  to  multiple  DMs.  each  with  its  own  LQ  objective  functional,  leads  to  the  feedback  Nash 


equilibrium  concept  [7],  Here,  one  finds  that  the  feedback  Nash  equilibrium  solution,  if  it  exists,  is 


obtained  by  solving  a  set  of  coupled  Riccati  equations.  However,  the  numerical  solution  of  these 


coupled  Riccati  equations  is  substantially  more  difficult  to  characterize.  Furthermore,  conditions 


insuring  existence  and  uniqueness  of  fixed  points  of  these  coupled  equations  have  not  been 


produced.  Nevertheless,  several  results  concerning  the  general  LQ  continuous-time  problem  have 


appeared  [7-12].  By  contrast,  the  discrete-time  case  has  received  considerably  less  attention. 


Indeed,  most  of  the  results  on  discrete-time  feedback  Nash  solutions  may  be  found  in  [7]  or  [13] 


But  these  authors  point  out  the  need  for  research  into  the  computational  aspects  of  solving  these 


game  problems. 


In  this  thesis,  we  solve  for  the  so-called  linear,  state-feedback  (perfect-state  information) 


Nash  equilibrium  of  a  discrete-time  descriptor  system.  For  simplicity  only  the  two  DM  case  is 


considered.  The  generalization  to  three  or  more  DMs  remains  an  open  problem.  The  material 


presented  here  suggests  an  obvious  approach  towards  extending  the  theory.  The  main  contribution 


of  this  work,  however,  is  the  development  of  algorithms  and  software  suitable  for  solving  coupled 


discrete-time  Riccati  equations  arising  from  LQ  feedback  Nash  games.  Additionally,  solutions  to 
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large-scale  problems  and  single-player  (control)  problems  follow  immediately  from  the 
descriptor-variable  formulation  taken  here.  The  convergence  behavior  of  these  equations  is 
investigated  which  results  in  a  contraction  mapping  argument  guaranteeing  existence  and 
uniqueness  of  a  solution  to  the  infinite-horizon  Nash  game  problem.  Also,  a  new  type  of  game 
theory  called  multirates  is  discovered  via  asymptotic  analysis.  The  tools  of  asymptotic  analysis  are 
also  used  in  producing  the  contraction  mapping  result.  The  descriptor-variable  game  formulation 
as  well  as  multirate  game  theory  have  not  been  examined  until  now.  These  theoretical 
contributions  enable  studies  of  a  significantly  larger  class  of  LQ  Nash  game  problems.  Moreover, 
there  are  several  practical  situations  where  this  new  theory  yields  more  accurate  and/or 
numerically  appealing  models. 

Although  theory  is  an  integral  part  of  this  work,  the  end  product  is  software.  Hence, 
software  engineering  is  a  key  issue.  That  is.  the  design  and  implementation  of  the  Riccati  software 
ought  to  comply  with  the  standards  anH  practices  currently  used  for  software  development.  This 
approach  ultimately  assures  the  quality  of  the  final  package.  However,  before  embarking  on  a 
detailed  discussion  of  the  algorithms  and  the  software  structure,  it  is  necessary  to  expand  on  this 
last  idea  more  fully. 


1.1.2  Computer-Aided  Control  System  Design 


A  study  such  as  this  falls  under  the  heading  of  Computer-Aided  Control  System  Design 
(CACSD).  Essentially,  this  research  field  strives  to  provide  the  control  system  and  related 
communities  with  high-quality,  reliable,  numerically  robust  algorithms  and  software.  CACSD  is 
still  a  relatively  young  area  of  research  (about  5  years  old).  This  remark  is  supported  by  the 
observation  that  formal  conferences  on  CACSD  are  relatively  new  [e.g..  14-16].  There  are  several 
problems,  such  as  linear  least  squares  or  generalized  eigenvalues-eigenvectors.  that  provide  a  basis 
for  solving  more  complicated  control,  system,  and  estimation  problems.  It  is  extremely  important 
that  stable  numerical  methods  are  used  for  solving  these  simple  problems.  Then  solutions  to  more 
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complicated  problems  rest  on  a  firm  algorithmic  foundation.  Some  tools  of  modem  numerical 
analysis  which  are  finding  significant  utility  in  CACSD  include  orthogonal  transformations  [17], 
Householder  transformations  [l7j.  Singular  Value  Decomposition  (SVD)  [17.18],  invariant 
imbedding  techniques  [19],  and  descriptor  variable  formulations  [19.20],  Some  of  these  tools  are 
useful  for  obtaining  matrix  forms  with  special  structure  (e.g..  block  triangular,  Hessenberg.  real 
Schur  form)  while  others  circumvent  difficulties  encountered  in  poorly  or  ill-conditioned  problems. 
As  a  consequence  of  the  CACSD  approach,  one  can  arrive  at  the  solution  to  a  problem  in  a 
straightforward  and  computationally  efficient  manner  by  exploiting  the  structure  imposed  by  a 
particular  technique.  A  recent  and  comprehensive  survey  of  the  preceding  ideas  may  be  found  in 
[21]. 

With  regards  to  robust  numerical  software,  it  is  widely  recognized  [e.g..  21]  that  the 
E1SPACK  [22.23]  and  UNPACK  [24]  software  packages  are  well-suited  for  generalized  eigenvalue- 
eigenvector  problems  and  linear  equation  problems,  respectively.  Both  packages  are  coded  in  the 
FORTRAN  [25]  programming  language.  Furthermore,  each  package  has  demonstrated  numerical 
superiority  over  the  years.  Hence,  this  software  is  a  natural  starting  point  for  building  the  coupled 
Riccati  equation  algorithms.  Sometimes  a  CACSD  package  consists  of  a  collection  of  subroutines, 
often  using  ElSPACk  and/or  UNPACK  as  the  lowest  level  routines.  RICPACK  [26],  a  software 
package  for  single  algebraic  Riccati  equations,  is  an  example  of  a  CACSD  package  with  this 
structure.  Other  times  a  more  substantial  undertaking  yields  a  package  capable  of  solving  a  broad 
range  of  problems. 

Despite  the  current  evolutionary  state  of  affairs,  several  trends  are  apparent.  For  instance, 
the  large  CACSD  software  packages  emerging  today  mav  be  classified  according  to  the  following 
groups:  menu-driven,  command-driven,  expert  systems,  anti  languages  Although  these  categories 
are  distinct,  examples  of  packages  can  be  fount!  that  possess  elements  ol  more  than  one  group. 
Because  of  the  volatility  of  CACSD  packages,  specific  examples  of  each  type  are  difficult  to  produce 
without  dating  this  text.  Nevertheless,  a  good  overview  ol  sollware  packages  available  today  may 
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be  found  in  [27],  In  particular,  we  mention  that  SIMNON  [28]  is  representative  of  a  command- 
driven  package  whereas  L-A-S  [27.  pp.  243-261]  qualifies  as  a  CACSD  language. 

Each  CACSD  software  group  has  its  strengths  and  weaknesses  in  terms  of  time  invested  by 
the  user.  For  example,  menu-driven  packages  have  the  advantage  that  the  infrequent  user  will 
probably  spend  a  minimal  amount  of  time  (re)learning  how  to  interact  with  the  package  by  virtue 
of  the  menu-driven  environment.  However,  the  disadvantages  include  the  fact  that  working 
through  pages  of  menus  is  ultimately  time-consuming.  More  importantly,  if  a  solution  procedure 
does  not  exist  as  one  of  the  choices  on  the  menus,  then  the  problem  is  quite  likely  unsolvable  by  the 
given  package.  On  the  other  extreme.  CACSD  languages  typically  require  much  more  time  (e.g.. 
hours)  to  (re)gain  familiarity  with  the  package.  But.  once  mastered,  an  almost  limitless  class  of 
problems  may  be  studied  depending  on  the  richness  of  the  language. 

Motivated  by  the  desire  to  conduct  systematic  numerical  studies  of  LQ  Nash  games  in 
discrete-time  as  well  as  the  need  to  efficiently  manipulate  matrices  in  a  user-friendly  environment, 
the  decision  is  made  to  integrate  the  coupled  Riccati  software  into  the  L-A-S  language  [27.  pp. 
243-261].  In  order  to  explain  the  subtleties  of  the  software  implementation,  a  brief  review  of  the 
L-A-S  language  package  is  required. 

1.1.3  The  L— A— S  Language 

BASIC.  FORTRAN,  and  PASCAL  are  standard  programming  languages.  Each  possesses 
qualities  and  attributes  that  are  characteristic  of  almost  any  ordinary  programming  language  in  use 
today  (e.g..  subroutine  capabilities).  It  is  desirable  that  a  CACSD  language  parallels  the 
organizational  model  set  by  these  familiar  and  well-established  computer  languages.  Furthermore, 
for  control  and  linear  system  problems  a  sophisticated  matrix  environment  is  mandatory.  In 
addition,  frequency  domain  techniques  require  the  analysis  and  manipulation  of  matrices  of 
polynomials.  It  is  according  to  these  prerequisites  that  the  L-A-S  language  was  created. 


L-A-S  stands  for  Linear  Algebra  and  Systems.  Furthermore.  L-A-S  is  a  CACSD  language  in 
the  strict  sense  of  the  word.  That  is,  L-A-S  conforms  to  standard  computer  science  definitions  for 
the  syntactical  specification  of  a  programming  language  [see.  27.  pp.  243-261.  and  29].  In  addition. 
L-A-S  has  been  tested  by  industry  and  academia  for  over  a  decade  at  over  a  dozen  locations  around 
the  world.  Numerically  speaking.  L-A-S  is  based  upon  the  EISPACK  [23]  and  LINPACK  [24] 
software.  Also,  the  NCAR  [30]  graphics  package  is  employed  to  provide  2D  and  3D  plotting 
capabilities.  In  summary,  there  is  ample  evidence  [27.29,31-36]  available  to  support  the  claim  that 
L-A-S  is  a  bona  fide  CACSD  language. 

In  the  normal  interactive  mode,  the  user  types  statements  directly  in  the  L-A-S  language 
interpreter.  Each  statement  is  either  a  command  to  the  interpreter  (e.g..  put  L-A-S  into  program 
debug  mode)  or  a  request  to  perform  some  kind  of  calculation.  The  former  instructs  L-A-S  to 
display  or  modify  various  status  information  concerning  the  current  L-A-S  work  session.  The 
latter  invokes  the  L-A-S  language  parser  which  subsequently  calls  upon  the  FORTRAN  subroutines 
needed  to  process  the  desired  computation. 

The  fundamental  concept  behind  any  L-A-S  statement  is  the  L-A-S  operator .  Essentially, 
operators  combine  input  data,  perform  some  desired  calculation,  and  generate  output  data.  The 
utility  of  L-A-S  operators  as  algorithmic  "building  blocks"  has  been  established  [27.35.36].  Thus, 
even  though  a  single  operator  may  not  be  available  to  solve  a  particular  problem,  it  is  quite  likely 
that  the  desired  result  can  be  obtained  by  concatenating  several  "lower-level"  operators.  The  L-A-S 
operators  are  divided  into  five  groups:  Input/Oulput.  Data  Handling,  Linear  Algebra.  Control 
Systems,  and  L-A-S  Program  Control.  Presently,  there  are  more  than  100  L-A-S  operators.  Also, 
the  user  may  define  up  to  100  matrices  with  the  total  number  of  matrix  elements  not  exceeding 
50.000.  The  maximum  order  of  any  particular  matrix  is  not  explicitly  limited. 

L-A-S  programs  are  written  by  combining  one  or  more  operator  statements.  Should  questions 
arise,  an  extensive  on-line  help  facility  containing  detailed  information  about  L-A-S  language  usage 
is  at  the  disposal  of  the  user.  In  totality.  L-A-S  and  its  supporting  software  consist  of  over  20.000 


7 


lines  of  FORTRAN  code  C 1977  Standard). 


The  L-A-S  language  will  be  used  extensively  throughout  this  dissertation.  Therefore,  it  is 


assumed  that  the  reader  is  adequately  familiar  with  L-A-S  so  that  the  L-A-S  programs  presented 


here  are  easily  understood. 


This  section  discusses  the  actual  computational  problem  studied  in  this  thesis  and  the  specifics 


of  the  contribution  of  this  work.  Because  the  coupled  discrete-time  Riccati  equations  analyzed  here 


are  deeply  rooted  in  Nash  game  theory,  two  game  scenarios  which  lead  to  the  solution  of  coupled 


discrete-time  Riccati  equations  are  developed.  First,  an  exposition  on  descriptor-variable  Nash 


games  is  presented.  Then,  multirate  descriptor  Nash  games  are  introduced.  The  mathematical  rigor 


associated  with  each  problem  is  deferred  until  Chapter  2.  The  purpose  of  this  discussion  is  to 


elucidate  the  theoretical  novelty  as  well  as  the  practical  applicability  of  descriptor  games.  In 


particular,  multirate  games  are  extremely  useful  for  formulating  optimization  problems  involving 


digital  communication  channels  operating  at  different  rates. 


Next,  the  details  of  the  thesis  contribution  are  highlighted.  The  computational  obstacles 


pertaining  to  iterating  two  coupled  discrete-time  Riccati  equations  are  delineated.  The  procedure 


by  which  they  are  overcome  is  outlined  and  justified.  Relevant  theoretical  issues  (such  as 


convergence  in  the  limit  as  the  number  of  iterations  tends  toward  infinity)  and  numerical  issues 


(such  as  preserving  symmetry)  are  addressed.  Finally,  the  software  implementation  is  described. 


Since  the  Riccati  software  is  integrated  into  the  L-A-S  language,  additional  care  must  be  taken  to 


insure  that  the  top-level  Riccati  routines  conform  to  the  L-A-S  interfacing  protocol.  Also,  the 


low-level  routines  must  be  engineered  properly.  Hence,  structured  programming  and  modularity 


concepts  are  discussed. 
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In  this  subsection,  we  describe  two  distinct,  yet  related  problems  where  coupled  discrete-time 


Riccati  equations  arise.  Both  problem  formulations  represent  new  contributions  to  the  theory  of 


Nash  games  utilizing  perfect  state  information.  However  the  fact  that  both  problems  lead  to  the 


solution  of  two  coupled  Riccati  equations  is  the  main  reason  for  their  inclusion  in  this  thesis.  The 


ensuing  discussion  is  primarily  qualitative. 


1 .2.1.1  Descriptor  Games 


A  linear  shift-invariant  (LSI)  discrete-time,  descriptor  system  takes  the  form  : 

E  x  (k+1)  =  A  x  (k)  +  B  i  u  ^k)  +  B2  u;(k)  (k  ^  0,  E  x  (0)  =  £  x0)  (1.2.1) 


yt(k)  =  Ci  x(k) 


(1.2.2a) 


y2(k)  =  C2  x  (k)  • 


(1.2.2b) 


This  game  problem  has  two  decision  makers.  DM1  and  DM2.  The  discrete-time  dynamic 


system  is  evolving  at  a  rate  indexed  by  the  integer-valued  variable,  k.  x(k)  €  R*  .  u^k)  €  RPt  . 


u2(k)  €  Rri .  yi(k)  6  Rm 1  ,  and  y2(k)  €  R1"2 .  The  matrix  E  is  square  and  is  assumed  to  be 


nonsingular  to  numerical  precision.  The  standard  state-space  formulation  is  recovered  by 


multiplying  (1.2.1)  by  E~^ .  The  system  being  described  is  depicted  in  Figure  1. 


Actually,  the  reasons  for  choosing  a  descriptor-variable  system  are  more  compelling  than  is 


first  apparent.  To  begin  with,  many  physical  systems  undergo  a  modelling  phase,  during  which 


lime  the  physical  laws  of  nature  are  applied  to  the  problem.  Often,  the  result  of  this  process  is  a 


static  and/or  dynamic  (i.e.,  algebraic  and/or  difference  equations)  descriptor-variable  description  of 


the  system  plant.  The  matrix  E  plays  the  role  of  a  mass  matrix  (when  Newton's  F  =  MA  is 


applied),  or  a  sparse  interconnection  matrix  (for  distributed  parameter  systems),  or  a  very  singular 


matrix  (for  economic  dynamic  games).  Whatever  the  case,  it  is  not  desirable,  or  even  feasible,  to  go 


to  the  standard  state-space  description  by  inverting  E.  Second,  the  theoretical  aspects  of 


descriptor-variable  systems  are  a  relevant  issue.  Allowing  E  to  be  singular  ultimately  enlarges  the 
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«,( k) 


u2(k) 

Figure  1.  Two-Player.  Discrete-Time,  Descriptor  Game. 


class  of  problems  that  can  be  studied  by  game  theorists.  Last  but  not  least,  the  descriptor-variable 
formulation  is  numerically  superior  to  state-space  ones  for  two  reasons.  Inverting  E  may  be 
numerically  impossible  or  inversion  of  E  would  destroy  any  inherent  structure  (e.g..  sparsity)  in 
that  matrix.  These  facts  and  more  support  the  philosophy  that  descriptor-variable  formulations  of 
optimal-control  and  dynamic-game  problems  are  physically,  theoretically,  and  numerically 
superior  to  state-space  formulations. 

Each  DM  has  an  associated  LQ  cost  functional  which  is  to  be  minimized.  We  solve  for  the 
linear,  state-feedback  (perfect-state  information)  Nash  equilibrium  solution  which,  by  definition, 
obeys  the  principle  of  optimality.  This  involves  extending  the  well-known  single-player 
optimization  result  to  descriptor  games.  This,  in  turn,  leads  to  the  solution  of  two  coupled 
discrete-time  Riccati  equations. 


1.2. 1.2  Multirate  Descriptor  Games 


This  subsection  introduces  the  theory  of  multirate  descriptor  games.  The  idea  is  that  one  DM 
is  constrained  to  control  the  system  at  a  rate  which  is  slower  than  the  other  DM.  It  turns  out  that 
multirate  descriptor  games  can  be  formulated  and  solved  via  single-rate  descriptor  game  theory  if 
appropriate  limiting  arguments  are  constructed.  In  addition,  there  are  several  practical  situations 
where  the  optimization  problem  is  more  accurately  modelled  by  multirates  than  by  single  rates. 

Singular  perturbation  techniques  are  used  successfully  to  exploit  the  presence  of  time  scales 
within  continuous-time  [37]  and  discrete-time  systems  [38],  These  applications  deal  with  time 
scales  that  are  inherent  in  the  underlying  system  plant  (i.e..  the  system's  eigenvalues).  By  contrast, 
relatively  little  attention  has  been  paid  to  the  case  where  time  scales  are  introduced  by  the  control. 
Consider  such  a  case  in  a  multiple-decision,  discrete-time  setting  where  each  DM  is  constrained  to 
control  the  system  at  a  different  rate.  Although  an  analysis  tool  has  been  developed  for  discrete¬ 
time  systems  with  multirate  samplers  [6],  it  is  not  suitable  for  multiple-decision,  optimization 
problems. 

Multirate  descriptor  games  begin  with  the  LSI  discrete-time,  noncooperative  game  problem  of 
the  last  subsection  and  constrain  the  D.Vls  to  play  at  different  rates.  As  before,  we  determine  the 
feedback  (perfect  state)  Nash  solution  [7]  which  leads  to  the  periodic  solution  of  two  coupled 
discrete-time  Riccati  equations.  Again,  we  consider  only  two  DMs.  However,  the  results  can  be 
extended  to  multiple  DM  situations  as  well.  In  addition,  we  restrict  attention  to  those  rates  which 
are  related  to  all  others  by  a  positive  integer  constant.  It  is  straightforward  to  generalize  to  rates 
which  are  related  by  rational  constants.  The  fast  player  (DM1)  plays  at  a  rate  that  is  an  integral 
multiple  of  the  rale  of  the  slow  player  (DM2).  Also,  the  control  policy  of  the  slow  player  relative 
to  the  fast  player  is  assumed  to  be  all-digital  in  that  the  slow  player  applies  zero  or  no  control 
during  those  instants  when  the  fast  player  is  acting  on  the  system  alone. 

Such  situations  arise  in  practice.  For  example,  consider  a  decentralized  control  problem  where 
the  controllers  must  communicate  with  the  system  via  a  digital  channel  and  the  maximum 


throughput  rate  of  each  channel  is  different.  A  second  application  of  multirate  descriptor  games 
occurs  in  the  field  of  agricultural  economics.  Here,  the  suppliers  typically  act  upon  the  market  at  a 
slow  rate  (e.g..  annually),  whereas  the  consumers  act  upon  the  market  at  a  fast  rate  (e.g..  daily). 
There  are  other  examples  of  problems  that  are  more  accurately  modelled  with  multirates  rather 
than  with  single  rates.  Moreover,  there  are  other  interesting  features  of  multirate  games  that  are. 
in  general,  not  present  in  single-rate  discrete-time  games.  For  a  complete,  self-contained  treatment 
of  multirate  Nash  game  theory  using  a  state-space  formulation,  the  reader  is  referred  to  [39]. 

1.2.2  Computational  Issues 

The  numerical  aspects  of  iterating  coupled  discrete-time  Riccati  equations  are  the  topic  of  this 
thesis.  The  computational  details  are  highly  nontrivial  because  of  the  coupling.  This  fact  is  the 
source  of  all  numerical  hardships  encountered  in  this  problem.  In  order  to  produce  any  algorithm 
for  solving  these  equations,  the  coupling  must  be  removed. 

Initially  the  relevant  equations  are  gathered  together  and  preliminary  notation  is  defined.  A 
cross-substitution  procedure  begins  the  decoupling  process.  Eventually  the  equations  are  rewritten 
in  a  form  that  permits  iteration.  At  this  point,  an  algorithm  is  presented  which  simultaneously 
iterates  the  Riccati  equations  and  solves  for  the  feedback  matrices  needed  if  the  problem  is 
motivated  by  a  descriptor  game.  As  given,  the  algorithm  is  designed  to  solve  single-rate  and 
multirate  descriptor  games  via  dynamic  programming.  But  this  fact  notwithstanding,  the  task  of 
iterating  a  set  of  coupled  discrete-time  Riccati  equations  is  still  accomplished.  Furthermore,  the 
conditions  which  govern  the  existence  of  a  solution  to  the  iteration  problem  (i.e..  existence  of  the 
next  Riccati  iterates)  become  apparent  and  are  stated  formally. 

However,  it  is  necessary  to  study  the  computational  aspects  of  iterating  coupled  Riccati 
equations  in  greater  detail.  Since  these  equations  involve  several  positive-(semi)definile.  symmetric 
matrices  and  several  quadratic  forms,  it  is  most  wise  to  seek  an  expression  where  these  quantities 
appear  explicitly  and  often.  With  this  goal  in  mind,  the  Riccati  equations  are  rewritten  in  a  form 
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more  amenable  to  computer  implementation.  Then,  the  specifics  of  the  coupled  Riccati  algorithm 
are  investigated.  Key  matrices  (usually  positive-(semi)definite  and  symmetric)  are  identified  and 
an  algorithm  for  computing  each  one  is  presented.  The  UNPACK  software  is  chosen  for  the  task  of 
manipulating  (i.e.,  decomposing,  inverting,  et  cetera)  these  key  matrices  while  exploiting  their 
special  structure  whenever  possible.  Not  surprisingly,  a  reduction  in  the  amount  of  computation 
results  from  this  method.  These  low-level  algorithms  are  subsequently  used  to  build  the  coupled 
Riccati  package.  This  structured  programming  technique  yields  highly  modular  and  efficient  code. 
These  aspects  of  the  software  engineering  process  make  the  coupled  Riccati  software  developed  here 
superior  to  other  approaches. 

Finally,  existence  and  convergence  issues  are  addressed.  It  is  natural  to  ask  if  the  coupled 
Riccati  iterates  converge  to  a  fixed  point  in  the  limit  as  the  number  of  iterations  tends  toward 
infinity.  To  date,  neither  necessary  nor  sufficient  conditions  insuring  such  convergence  have  been 
established.  This  work  investigates  existence  issues  associated  with  finite  horizon  problems  and 
convergence  issues  associated  with  infinite  horizon  problems.  In  particular,  a  contraction  mapping 
argument  is  developed  that  guarantees  existence  of  and  convergence  to  a  unique  fixed  point  for  the 
infinite  horizon  coupled  Riccati  problem. 

The  main  analysis  tool  used  to  obtain  the  convergence  results  is  asymptotic  analysis.  First,  a 
small  parameter,  is  introduced  into  one  DM’s  cost  functional.  As  £i  — *  0.  the  two-player  LQ 
descriptor  Nash  game  problem  reduces  to  the  LQ  descriptor  regulator  problem.  Thus,  an  initial 
"bridge"  between  the  fields  of  optimal  control  and  Nash  games  is  established.  Furthermore,  by 
setting  €[  =  0  periodically,  a  new  game  called  multirates  is  created.  This  is  the  essence  of  the 
limiting  argument  mentioned  in  the  last  subsection.  Moreover,  there  are  several  practical  situations 
where  multirate  game  theory  is  more  applicable  than  standard  single-rate  game  theory. 

Next,  a  second  small  parameter.  e2.  is  introduced  into  the  other  DM  s  cost  functional.  Then 
we  let  €j  -» 0  and  e2  -*  0  independently.  Subsequently,  we  discover  that  under  appropriate 
assumptions,  there  exists  a  region  where  the  coupled  discrete-time  Riccati  iterations  behave  as  a 
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contraction  mapping.  Therefore,  we  produce  sufficient  conditions  that  guarantee  that  the  Riccati 
iterates  will  converge  to  a  unique  fixed  point  in  the  limit.  Furthermore,  empirical  simulations  have 
verified  the  contraction  mapping  behavior  and  indicated  that  the  bounds  obtained  are  rather 
conservative. 

1.2.3  Software  Implementation 

The  computer  algorithms  and  software  developed  for  solving  both  single  and  coupled 
discrete-time  Riccati  equations  are  integrated  into  the  L-A-S  CACSD  language  as  new  operators. 
This  approach  to  software  implementation  is  novel  and  unique.  Theory  [19.40]  and  software  [26] 
have  been  developed  for  the  numerical  solution  of  single  algebraic  Riccati  equations.  Additional 
work  that  is  very  closely  related  to  this  problem  includes  [41-44],  However,  until  now  similar 
efforts  for  coupled,  discrete-time  Riccati  difference  equations  have  not  been  undertaken. 
Furthermore,  descriptor  variable  theory  [19.20.43]  is  applied  to  the  LQ  Mash  game  formulation 
which  leads  to  the  so-called  generalized  coupled  Riccati  equations.  Thus,  the  class  of  problems  that 
can  be  solved  is  broadened. 

Each  new  operator  involves  careful  structuring  of  the  algorithm  as  well  as  a  sound  software 
engineering  basis  for  writing  the  codes.  The  high-level  Riccati  routines  interface  with  the  L-A-S 
protocol.  The  low-level  routines  are  highly  modular  and  computationally  efficient.  In  the 
computational  analysis  it  is  shown  that  the  task  of  iterating  coupled  Riccati  equations  reduces  to 
solving  systems  of  linear  algebraic  equations  where  one  or  more  matrices  have  special  properties 
(like  positive-definite  and  symmetric).  Hence,  the  decision  to  build  the  low-level  routines  from 
calls  to  the  UNPACK  library  is  fairly  justified.  Further,  the  high-level  routines  are  built  from 
frequent  calls  to  the  low-level  routines  as  good  structured  programming  practice  dictates.  The 
overall  software  structure  is  depicted  in  Figure  2. 

Altogether,  there  are  a  total  of  six  new  L-A-S  operators.  Their  mnemonic  names  are  SYST. 
I.Q.  DRE.  GAME.  LQNG.  and  V1LTR.  Collectively  they  form  a  single  and  coupled  discrete-time 
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Figure  2.  Coupled  Riccati  Software  Structure  and  Interface  with  L-A-S 


Riccati  equation  solver  subpackage  of  L-A-S.  Descriptor-variable  systems  are  handled  directly. 
Also,  auxiliary  codes  are  required  for  repeatedly  performing  small  tasks  (e.g..  multiplying  two 
general  matrices).  The  new  L-A-S  operators  plus  the  associated  auxiliary  support  routines  are 
written  entirely  in  FORTRAN  77  [25]  and  amount  to  approximately  3000  lines  of  code. 
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1.3  Thesis  Organization 

This  dissertation  is  organized  in  the  following  manner  Chapter  2  formulates  two  distinct 
yet  related  two-player.  LQ  dynamic  games  using  descriptor- variable  theory  The  first  has  DMs 
that  control  the  descriptor  system  at  a  single  rate  while  the  second  restricts  one  DM  to  apply 
control  at  a  slower  rate.  Assuming  perfect-state  information,  linear  state-feedback  Nash  strategies 
are  determined  for  each  case.  The  solution  to  both  problems  requires  the  iteration  of  two  coupled 
discrete-time  Riccati  equations.  The  purpose  of  this  chapter  is  three-fold.  First,  bas.  cepts  and 
notation  are  introduced.  Second,  the  coupled  Recall  problem  is  motivated.  Third,  the  theory  of 
descriptor  games  and  multirate  games  is  formalized  and  documented. 

Chapter  3  investigates  the  computational  aspects  associated  with  iterating  two  coupled 
discrete-time  Riccati  equations.  A  cross-substitution  procedure  is  used  to  decouple  the  equations 
A  preliminary  iteration  algorithm  is  presented  in  a  game  context.  Then  attention  is  briefly  directed 
to  existence  issues  (iterability  -  the  ability  to  iterate  a  recursive  equation)  Next  the  coupled 
Riccati  equations  are  rewritten  in  a  form  more  amenable  to  computer  implementation  Matrices  are 
identified  that  possess  special  structure  (e  g.,  positive/symmeinc)  and  recur  repeatedly  throughout 
the  equations.  Algorithms  for  computing  these  matrices  are  given.  Next,  calculations  of  the 
feedbacks  and  the  Riccati  iterates  are  described.  Then,  the  multirate  Riccati  algorithm  is  presented 

The  existence  of  solutions  to  finite-horizon  problems  and  convergence  of  Riccati  iterates  for 
infinite-horizon  problems  are  studied  in  greater  depth.  Several  new  results  are  stated.  The 
contraction  mapping  argument  is  developed  here. 

Chapter  4  is  devoted  to  the  new  L-A-S  operators  created  for  solving  single  and  coupled  Riccati 
iterations.  Details  of  the  design  and  syntax  of  each  operator  are  given.  Also,  examples  illustrate 
how  the  L-A-S  language  may  be  used  to  study  single-rate  and  mullirate  LQ  descriptor  Nash  games. 
Chapter  5  summarizes  this  work  and  discusses  future  research  topics. 

Three  appendices  are  included  which  support  the  theoretical  and  numerical  results  presented 
in  the  main  body  of  this  thesis.  Appendix  A  contains  selected  software  listings  of  the  low-level 
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Riccati  routines.  Appendix  B  consists  of  facts  and  lemmas  that  are  used  in  proving  the  contraction 
mapping  result.  Appendix  C  contains  the  L-A-S  program  run  which  produced  the  data  used  in  the 
contraction  mapping  example  of  Subsection  3.4.3. 


I 

£ 
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CHAPTER  2 
PRELIMINARIES 


This  chapter  formulates  and  solves  two  LQ  dynamic  game  problems  in  discrete  time.  The 
first  problem  considers  a  descriptor-variable  system  where  the  DMs  control  the  plant  at  the  same 
rate.  The  second  problem,  a  multirate  game,  further  restricts  one  DM  to  apply  control  at  a  slower 
rate.  Assuming  perfect-state  information,  linear  state-feedback  Nash  strategies  are  determined  for 
each  situation.  Of  more  interest,  however,  is  the  fact  that  both  problems  lead  to  the  solution  of 
two  coupled  discrete-time  Riccati  equations. 

Although  the  descriptor-variable  formulation  has  been  applied  to  the  LQ  regulator  problem 
[2 1  ].  it  has  not  been  attempted  for  dynamic  games.  Hence,  this  approach  is  completely  new. 
Moreover,  multirate  descriptor  game  problems  are  heretofore  posed  yet  unsolved.  As  such,  this 
chapter  additionally  provides  significant  extensions  to  the  theory  of  infinite  dynamic  games. 

However,  the  primary  purpose  of  the  following  discussion  is  to  motivate  the  coupled  Riccati 
problem  whose  computational  aspects  are  studied  in  the  remainder  of  this  thesis.  Second,  basic 
concepts  and  notation  are  defined. 

2. 1  LQ  Discrete— Time  Descriptor  Nash  Games 

Subsection  2  1.1  formulates  an  LQ  descriptor  game  in  discrete  time.  Subsection  2.1.2  states 
the  Nash  equilibrium  solution  to  the  proposed  game. 


2.1.1  Problem  Formulation 


%• 

*»  11 

Consider  the  linear  shift-invariant  discrete-time,  descriptor  system  described  by 

< 

E  .t  (k+1 ) 

=  A  -V  ( k )  +  Bx  u  i  ( k )  +  #2u2(k) 

( k  ^  0.  E  x  (0)  =  E 

(2.1.1) 

k 

v ,( k ) 

=  C ,  x  ( k ) 

(2.1.2a) 

s 

y  2(  k ) 

=  C  2  x  ( k )  • 

(2.1.2b) 

■ 


18 


This  game  problem  has  two  decision  makers,  DM1  and  DM2.  The  discrete-time  dynamic 
system  is  evolving  at  a  rate  indexed  by  the  integer-valued  variable,  k.  x  (k)  €  Rm  .  Ui(k)  €  Rfl , 
u2(k)  €  Rfi  ,  y/k)  €  Rmt  .  and  y2(k)  €  Rm* .  The  matrix  E  is  square  and  is  assumed  to  be 
nonsingular  to  numerical  precision.  Multiplying  (2.1.1)  by  E~ *  yields  the  standard  state-space, 
dynamic  game  formulation. 

Each  DM  must  have  an  objective,  which  may  or  may  not  agree  with  the  other  DM's  objective. 
Assume  that  an  LQ  cost  functional  has  been  selected  for  assessing  the  payoffs/losses  incurred  by 
each  DM.  This  is  frequently  done  in  modern  control/game  applications.  Then,  the  performance 

criterion  to  be  minimized  by  each  DM  is  given  by 

A  ,  r, 

J  lCyi-V’)  =  y  y[(k+l)  S/k+l)  y  i(k  +  l)  +  u[  (k)  R  x(k )  u  ,(k)  (2.1.3a) 

^  k=o 

J3(y i.y2)  t  y  £.  >'2  (k+1)  S2(k+1)  v2(k+l)  +  (k)  ^2(k)  u2(k)  (2.1.3b) 

2  lc=0 

where  y2  denote  the  mappings  from  the  information  set  to  the  control  and  T  denotes 
transposition.  The  time-varying  matrices  S,(k+l)  and  l?,(k).  i=1.2  are  symmetric  and  positive 
semidefinite.  Equation  (2.1.3)  describes  a  finite-time  game  of  duration  Tf  .  If  Tf  -*  00,  then  the 
infinite-lime  problem  can  be  studied.  Let 

H  J  (0.  1.  2 . Tf  I  (2.1.4) 

denote  the  index  set  of  the  horizon  of  the  problem.  Also,  let  N  ^  { 1 .  2)  be  the  set  used  to  index  the 

decision  makers.  Next,  recall  the  definition  of  a  Nash  equilibrium  for  a  two-player  game  : 
Definition  2.1  :  Nash  Equilibrium  [see  7  for  details] 

A  set  of  strategies  f  y[,  y2]  constitutes  a  Nash  equilibrium  solution  if  and  only  if 

/,*  t  J  i(y,*.y2’)  ^  J  ,(y , .  y2’)  (2.1.5a) 

J2‘  t  72(y,’.  y2*)  <  y2(yj*.  y2)  (2.1.5b) 

for  all  {  y^  y2l  €  Cl  -  the  set  of  all  admissible  mappings.  A  star  denotes  the  value  of  that  quantity 

at  the  Nash  equilibrium  solution. 


In  this  chapter,  we  solve  for  the  so-called  linear,  state-feedback  (perfect-state  information) 
Nash  equilibrium  [7].  By  definition,  this  solution  obeys  the  principle  of  optimality  Since  each  DM 
will  assume  that  the  other  DM  will  play  a  linear,  state-feedback  strategy,  the  feedback  (perfect- 


state)  Nash  policy  of  each  player  in  this  game  is  defined  by 


u  jCk)  : 

=  yMx(k)) 

=  — /’^k)  x  (k) 

(2.1.6a) 

u2(k)  : 

=  y2k  (x  (k)) 

=  — F2(k)x(k)  . 

(2.1.6b) 

Then. 

Vi  = 

Vi°.  Vi1 . 

(2.1.7a) 

>2  = 

y2  ■  V21 . 

yL  ]  • 

(2.1.7b) 

We  solve  the  basic  LQ  game  problem  described  above  in  the  next  subsection. 


2.1.2  Nash  Equilibrium  Solution  of  LQ  Discrete— Time  Descriptor  Games 


Define  the  system  matrix  of  each  DM  as  follows  : 


A  ](k)  t.  A  -B2F2( k)  and 


-4  2(  k )  “  A  -  S,F,(  k)  . 


(2.1.8a) 

(2.1.8b) 


Then  utilizing  (2.1.1 )-( 2.1.2)  and  (2.1.6)  as  well  as  (2.1.8).  we  can  write  (2. 1.3)  as 


T , 


j  i(Vi.‘y:)  =  Y  T.  *r  (k) 


k=0 

Tf 


\AcL  (k)  C  jS^k+1  )C ,  Acl  Ck )  +  F\ (k)  R  ,(k)  F,(k) 


x  (k)  (2.1.9a) 

x(k)  (2.1.9b) 


J2(y i.y2)  =  y  22  x '  ( k )  ^/£(k)  C2S2(k+l  )C2  ACi  (k)  +  F2  (k)  F2(k)  F2(k)  | 
where  .4,  L  (k)  ^  4  j(k)  —  B  ,F,(k)  =  .4  2(k)  —  B2F2(k).  Here  the  subscript  CL  means  closed 

loop. 


Because  the  feedback  (perfect  state)  Nash  equilibrium  solution  obeys  the  principle  of 
optimality,  this  two-player  game  is  equivalent  to  two  coupled  one-plaver  LQ  regulator  problems. 
The  following  result  concern. ,ig  single-player  optimization  has  been  extended  to  descriptor  systems. 
It  forms  the  basis  for  the  entire  coupled  Riccati  equation  study. 


j  Wi  *u  *  -  * 


Fact  2.1  :  Optimal  LQ  Regulator  for  Discrete-Time  Descriptor  Systems 

The  linear  state  feedback  policy  yv  i  €  N  defined  in  (2.1.6)  satisfies  (2.1.5)  subject  to  (2.1.1)- 


(2.1.3)  if  and  only  if 

u  j*(k)  =  —  (/? j(k)  +  BfJf^k+l)#,)-1  (£fJf,(k+lM,(k))  x*(k)  (2.1.10) 

where  fifi(k)  satisfies  the  following  generalized  discrete-time  Riccati  equation  : 

Et  K  fX)  E 

=  Af  (k)ifi(k+lMj(k)  —  (Af  (k)Xj(k+l)Bi)  (T^k))-1  (5ftf,(k+lM  ,(k))  +  Q,(k)  (2.1.11a) 

=  Af (k)  j  If  ,(k+l)  —  A'i(k+l)5i(ri(k))-15fA'i(k+l)  A,(k)  +  Q,(k)  (2.1.11b) 

=  >1 T (k)  [(A’i(k+l))“1  +  4,(11  t(k))r*£f  p1  >1  j(k)  +  £>i(k)  (2.1.11c) 


(2.1.11c) 


<2,(k)  Z  Cf  Sj(k+1)  Cj. 

r,(k)  t  RW  +  BfK^+DBi 

with  terminal  constraint 


Et  K  [T f  +1  )£  =  <2,(7.  )  =  Cf  5,(7,  +1)C,. 


Hence. 


F,(k)  =  (*,(k)  +  Bf/f^k+Dfi,)-1  (B[K  ,(k+lM  ,(k)) 
72(k)  =  (i?2(k)  +  S^2(k+1)52)*1  (5£A:2(k+lM2(k)) 


(2.1.12) 


(2.1.13) 


(2.1.14) 


(2.1.15a) 

(2.1.15b) 


Proof  :  For  the  case  E  —  I,  this  result  is  well  known.  A  proof  of  this  case  may  be  found  in  [7. 
p.22l].  The  extension  to  the  case  of  arbitrary  (and  possibly  singular)  E  has  been  investigated  in 
[43.44].  However,  it  is  illustrative  to  outline  the  major  components  of  the  proof.  Let  p,(k)  denote 
the  codescriptor  vector  for  the  i,h  player,  i  €  N.  Then  each  DM  faces  the  following  two-point 


boundary  value  problem 

E  00  x (k+1 ) 

0  Af  0  p,(k+l) 

0  -Bf  0  u  j(k+l ) 

with  boundarv  conditions 


A,  0  Bt  .r(k) 

~Q>Er  0  p,(k) 

0  0  R,  U,(k) 


(2.1.16) 


mm 


v;.'  ..v*;  wm 


iki 


E  x(0)  —  E  x0 


ET  P$Tf  )  =  ET  K-JLTf  )E  x(Tf). 

Hence,  the  relationship  between  descriptor  and  codescriptor  vectors  is 


(2.1.17) 


(2.1.18) 


/>i(k)  =  A'i(k)£  x(k)  .  (2.1.19) 

In  view  of  (2.1.5).  it  is  well-known  [7.45]  that  the  necessary  (and  sufficient)  conditions  for 

existence  of  a  minimum  of  (2.1.3)  restricted  to  the  class  of  linear  state  feedback  policies  relies  upon 

application  of  the  Matrix  Minimum  Principle  (MMP)  [46]  to  the  cost  functionals.  Noting  that,  in 

general,  the  elements  of  F  \  and  F 2  are  independent,  it  is  straightforward  to  take  the  state  feedback 

form  (2.1.9)  of  /j  and  J  2  in  order  to  calculate 

=  0  implies  u/(k)  =  -  U,(k)  +  iCk+Dfl,)-1  (BfKfk+DA  >(k))  x*(k)  .(2.1.20) 
Thus,  (2.1.10)  and  (2.1.15)  are  verified.  The  discrete-time  Riccati  equation  (2.1.11)  can  be  derived 
from  (2.1.16).  Let  G,  ^ ,_1  Bf.  Then  (2.1.16)  is  equivalent  to 


(2.1.21) 


£  G,  x (k+1 )  _ 

A;  0 

x  (k) 

0  Af  p,(  k+1) 

-QiET 

/>,(k) 

Therefore. 

£  x  (k+1)  +  G,  />j(k+l)  =  A  j  x  (k) . 
But  from  (2.1.19)  we  have 

I  A*,’1  +  G,  L,(k+1)  =  A  j  x  (k) . 


(2.1.22a) 


(2.1.22b) 


Af  p,(  k+1)  =  Af  K~l  +  G,  A  j  x  (k)  =  —Q,x(k)  +  ET  K,E  x(k).  (2.1.22c) 

Since  (2.1.22c)  must  be  true  for  all  x  (k).  we  require  that 

a  r  I  is -1  j.  n  1  *  a  j.  n  =  e-/ ir  v  a  i  'nn') 


Af  j  ATT1  +  G,  .4,  +  Q,  =  E1  A',£  .  (2.1.22d) 

It  is  clear  that  (2.1.11c)  and  (2.1.22d)  are  equal.  The  terminal  condition  (2.1.14)  comes  from 
setting  Tf  =  — 1  in  (2.1.3).  See  [5.7.44]  for  details. 
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Remarks  : 

1)  Equation  (2.1.11)  represents  several  forms  of  the  discrete-time  Riccati  equation.  All  are 
equivalent  assuming  the  appropriate  inverses  exist. 

2)  The  coupling  of  the  Riccati  equations  (2.1.11),  i  €  N  occurs  through  the  feedback  matrices. 
■Fj(k).  To  see  this,  substitute  (2.1.15)  into  (2.1.8).  The  resulting  equations  are  coupled.  Since 
(2.1.11)  depends  on  (2.1.8)  and  (2.1.15).  the  Riccati  equations  are  coupled  too. 

3)  After  the  lengthy  discussion  in  Chapter  1  about  the  advantages  of  descriptor-variable 
formulations,  it  is  fortunate  that  the  matrix  E~ *  appears  nowhere  in  the  derivation.  However, 
a  close  examination  of  (2.1.11)  reveals  that  A"j(k)  is  not  easily  obtainable  given  the  right-hand 
side  of  the  equation.  In  fact,  the  case  of  E  singular  poses  an  interesting  problem  because  then 
there  may  be  zero.  one.  or  multiple  Aj(k)'s  depending  upon  the  right-hand  side.  Nevertheless, 
the  Riccati  algorithm  presented  in  Chapter  3  will  recover  A,(k)  without  explicitly  inverting  E . 
Because  of  the  complications  which  arise  from  a  singular  E  matrix,  we  impose  the  following 
restriction. 

Assumption  2.1  :  Throughout  the  remainder  of  this  thesis,  the  matrix  E  is  always  assumed  to  be 
nonsingular  to  working  numerical  precision. 

2.2  Multirate  LQ  Descriptor  Nash  Games 

This  section  develops  the  theory  of  multirate  descriptor  Nash  games.  The  idea  is  that  one  DM 
is  constrained  to  control  the  system  at  a  rate  which  is  slower  than  the  other  DM.  It  turns  out  that 
multirate  games  can  be  formulated  and  solved  via  single-rate  Nash  game  theory  if  appropriate 
limiting  arguments  are  constructed.  Theorem  2.1  in  this  section  makes  the  last  statement  precise. 
In  addition,  there  are  several  practical  situations  where  the  optimization  problem  is  more 
accurately  modelled  by  multirates  than  by  single  rates.  Chapter  1  contains  those  details. 

The  organization  of  this  section  is  identical  to  the  last  one.  Subsection  2.2.1  formulates  the 
problem.  Subsection  2.2.2  discusses  the  Nash  solution  of  this  multirate  game  problem.  In  this 
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thesis,  the  multirate  theory  is  introduced  here  and  mentioned  in  subsequent  chapters  as  a  special 
topic.  The  reader  is  referred  to  [39]  for  a  complete,  self-contained  treatment  of  multirate  game 
theory  using  a  state-space  formulation. 


Consider  the  linear  shift-invariant  discrete-time  system  described  by  (2.1.1W2.1.2).  This 
game  problem  has  two  decision  makers.  DM1  and  DM2.  DM1  will  be  referred  to  as  the  "fast" 
player  and  DM2  will  be  referred  to  as  the  "slow"  player  throughout  the  discussion.  The  discrete¬ 
time  descriptor  system  is  evolving  at  a  rate  indexed  by  the  integer-valued  variable,  k.  DM1  plays 
at  rate  k.  i.e.,  at  every  k.  However  DM2  is  constrained  to  play  slower,  say  at  rate  j.  where  j  is  also 

an  integer-valued  variable  and  is  related  to  k  by  a  positive  integer  N  as  follows  : 

^  k/N  whenever  k  /  N  €{0.1.2....}  (2  2  1) 

^  ~  undefined  otherwise 

Thus,  j  simply  indexes  the  state  of  the  underlying  dynamic  system  (2.1.1  )-(2. 1.2)  sampled  every 
Nr/l  time.  But  j  and  k  are  really  dummy  variables.  Hence,  it  is  preferable  to  define  a  real-valued 
variable.  L  ^  k  /  N  .  for  all  k.  Then  joint  interaction  of  DM1  and  DM2  on  the  underlying  system 
occurs  if  and  only  if 


L  t  —  €  {0.  1.  2 _ }  t  Z+  .  (2.2.2) 

“  N 

The  condition  (2.2.2)  will  play  a  central  role  in  the  characterization  of  the  solution  to  this 
multirate  game  problem. 

The  performance  index  to  be  minimized  by  each  DM  is  given  by  (2.1.3)  as  before.  The  time- 
varying  matrices  Sj(k+1)  and  ^,(k).  i=l,2  are  symmetric  and  positive  semidefinite.  Equation 
(2.1.3)  describes  a  finite-time  game  of  duration  Tt  .  If  Tt  —  oo,  then  the  infinite  horizon  multirate 
LQ  descriptor  Nash  game  may  be  studied.  The  definitions  (2.1.6)  and  (2.1.7)  are  valid  as  well. 
Furthermore,  we  seek  the  feedback  Nash  equilibrium,  so  (2.1.8 )-( 2. 1.9)  hold.  However,  this  Nash 
equilibrium  will  be  for  a  multirate  game:  hence,  some  differences  in  problem  formulation  are  to  be 
expected. 
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In  order  to  account  for  the  fact  that  one  DM  is  applying  control  at  a  slower  rate,  we  define  the 
all-digital  control  policy  as  follows  : 


All  Digital  Control  Policy 


if  L  €  Z+ 

else  (no  control  applied) 


(2.2.3) 


Remark  :  In  some  applications,  a  multirate  Nash  game  may  arise  where  the  "slow"  player(s) 
employs  a  sampled-data  control  policy,  that  is.  the  control  law  is  held  at  its  previous  value  until  k 
is  such  that  the  condition  (2.2.2)  holds.  Such  problems  will  be  addressed  elsewhere. 

We  solve  the  basic  multirate.  LQ  game  problem  described  above  in  the  next  subsection.  In 
doing  so.  we  will  discuss  the  meaning  of  a  Nash  equilibrium  in  a  midtirate  setting. 

2.2.2  Nash  Equilibrium  Solution  of  Multirate  Descriptor  Games 

In  this  subsection,  we  show  that  under  appropriate  assumptions,  the  multirate  game  defined  in 
Subsection  2.2.1  can  be  solved  using  standard  single-rate  game  theory  with  minor  modifications. 
The  discussion  and  results  of  Subsection  2.1.2  are  valid  here.  Now.  we  construct  a  limiting 
argument  which  shows  that  as  the  small  parameter  €  -*  0,  the  all-digital  control  policy  (2.2.3)  is 
achieved  asymptotically. 

For  notational  reasons,  let  us  study  a  slightly  simplified  version  of  the  all-digital  multirate 
game  formulated  in  the  previous  subsection.  Consider  the  weighting  matrices  as  an  ordered  pair. 
Then,  we  are  interested  in  the  quantities  (S  i(k+l),  R  j(k))  and  (S2(k+1).  2(k)). 

Define  : 

(S  ](k  +  l ).  R  ]Ck ))  =  (Sp  R  [)  where  S ,  ^  0  and  Rj  >0  are  constant  matrices  for  all  k.  (2.2.4) 
and 

(S2,  R 2)  where  .V2  ^  0  and  R  2  >  0  when  L  €  Z+ 

(S2(k+1).  /?2(k))  *  •  (2.2.5) 

(o.  in  where  e-><)  otherwise 

€ 

Here.  0  is  the  zero  matrix  and  I  is  the  identity  matrix. 
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The  main  results  are  equally  applicable  to  arbitrary  but  bounded  S^k+l).  S2( k+1).  ^i(k). 
and  fl2(k),  (  when  L  6  Z+  ).  Furthermore,  this  setup  accounts  for  the  fact  that  DM2  is  applying 
control  at  a  slower  rate.  In  fact,  as  €-*0,  the  optimal  control  of  DM2  approaches,  in  the  limit,  the 
all-digital  control  policy  described  by  (2.2.3).  Furthermore,  under  mild  boundedness  assumptions, 
both  J  i  and  J2  are  finite.  Hence,  the  problem  is  well-posed.  These  last  statements  are  supported 
by  the  next  result. 

Theorem  2.1  :  Characterization  of  Multirate  Nash  Games 

For  the  pair  (S2(k+1).  i?2(k))  as  defined  in  (2.2.5).  if  /4j(k)  and  X^k+l).  1=1.2  remain  bounded 
for  all  k  €  H.  then  as  €-*0.  llu2(k)ll-*0  whenever  L  i  Z+.  Further,  whenever  L  £Z+  . 

II  A  J 1  (k)  II  — •  0  where 

A7,(k)  *  yf  (k+1  )S,(k+l)y,(k+l )  +  uf  (k)R,(k)u  ,(k)  „  (2.2.6) 

Moreover.  u2  (k)  isO(e)  and  so  is  A72  (k). 

Proof  :  Throughout  the  proof  all  norms  taken  are  assumed  to  be  the  standard  Euclidean  norm. 
Also,  it  is  useful  to  define  : 

<r(  )  ^  largest  singular  value  of  (  )  =  II.  and 
<r(  )  ^  smallest  singular  value  of  (  )  . 

It  is  given  in  the  problem  statement  that  B  j.  B2.  Si.  S2.  and  R  ]  are  matrices  with  bounded 

entries.  Now  as  €  — *  0.  II  £->(k)  II  =  —  —*  oo  whenever  L  t  Z+.  In  fact.  <j(R ->(k))  =  (FfAitk))  =  — 

e  ~  '  € 

by  construction.  For  the  moment,  assume  that  K  j(k  +  l ).  A'2(k+1 ).  A  j(k).  and  A  2(k)  exist  and  are 

finite.  This  requirement  is  investigated  further  in  the  next  chapter.  Consequently,  for  sufficiently 

small  €.  II  r2(k)  II  ^  »  oo.  Specifically,  from  the  relationship 

ll(X)~*l|  =  <f((X)-l  )  =  - ? -  for  anv  matrix  .Y  where  o;( -Y  )^0  .  (2.2.7) 

a(X) 

and  the  fact  that  lor  any  two  positive  semidefinite  symmetric  matrices  X  .  Y  . 


I  ( X  +Y  )-!  I  =  - - -  — - —  =  I  (X  )~l  I 

a(X+y)  a(X) 

we  conclude  that 


i  (r2(k)|  1 1  =  |s:(r2(k))  ]" 


1  * 1 

<  a(R2(k))  ~  =  I  =  €  -  0. 

e 


Hence. 


(2.2.8) 


I  u2  (k)  I  =  l-(*2(k)  +  BT2K2(.k+\)B2)-1  (5£X2(k+lM2(k))  x*(k)l 

=  l-(r2(k))"1  (5jX2(k+lM2(k))  x‘(k)  I 
^  l-(r2(k))-1l  •  I  fl  2X 2(k+ 1 M  2(k)  I  •  I  x*(k)  I 

=  e  I  B2K2(k+\)A  2(k)  I  •  lx*(k)l  -»  0  for  all  I  x*(k)  I  <  oo  . 


(2.2.9a) 

(2.2.9b) 

(2.2.9c) 


(2.2.9d) 


Thus.  u2(k)  is  O(e).  That  is.  there  exist  non-negative  constants.  M  j  and  M 2.  such  that 


A/,«  lu2*(k)l  <  M2e. 


example. 


choose 


M  \  —  0 


M2  =  K8^2(k+lM2(k))x'(k)I.  Obviously.  ^2(k)  and  r2(k)  are  O  —  . 


Having  established  that  u2  (k)  is  O  (e),  we  show  that  A J 2  (k)  is  also.  When  L  i.  Z+,  we  have 


from  (2.2.5)  and  (2.2.6)  that 


AJ2(k)  =  u £'(k)  £2(k)  u2  (k)  . 

Clearly.  A/2  (k)  is  O(e)  O III  0(e)  =  0(e). 


(2.2.10) 


Remark  :  From  Theorem  2.1.  we  conclude  that  lim  II  |r->(k)|  5  =  0  whenever  L  €  Z+  .  Hence 

e-n  t  *  I 

the  discrete-time  Riccati  equation  (2.1.11a)  reduces  to  the  discrete-time  Lyapunov  equation  : 


ET  K 2(k)  E  =  /1 2  (k)X2(k+l  M  2(k)  +  Q,( k) . 


Therefore,  as  €—*0,  II  u2  (k)  fi  as  well  as  1  A J2  (k)  II  -*  0.  and  the  all-digital  control  policy  of 
(2.2.3)  is  realized  in  the  limit.  Hence,  the  multirate  descriptor  Nash  game  problem  proposed  in 
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Subsection  2.2.1  can  be  completely  characterized  using  single-rate  theory,  if  the  pair 

(S2(k+l).  R 2(k))  =  (0,  —I)  and  e-*0  whenever  L  (Z+ .  Indeed,  the  multirate  Riccati  iterations 

6 

are  decoupled  when  L  i  Z+.  But  when  L  €  Z+.  the  coupling  is  present  as  in  the  first  problem  of 
Section  2.1.  Thus,  multirate  games  require  the  periodic  iteration  of  coupled  discrete-time  Riccati 
equations. 
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2.3  A  Numerical  Example 

In  this  section,  a  numerical  example  is  constructed  which  illustrates  that  there  are  situations 
where  a  descriptor-variable  game  formulation  is  computationally  superior  to  the  corresponding 
state-space  game  formulation  in  terms  of  numerical  accuracy.  Consider  the  system  (2.1.1  )-(2. 1.2) 
with  performance  indices  (2.1.3).  Let  C1=C2  =  /?i  =  #2  =  I-  Then  choose 


0.5  -1.0  -1.0 

9.8  7e  2 

1.37 

A  = 

0.0  0.77  -1.0 

.  B !  = 

1.23 

.  b2  = 

—l.Oe  -3 

0.0  0.0  0.001 

— l.Ole  -3 

l.Oe -5 

and 


E 


l.OeO  l.Oe  3  l.Oe  6 
0.0  9.9e  5  l.Oe  9 

— 1.95677e  — 1  0.0  l.Oe  15 


.Notice  the  wide  range  of  magnitudes  in  E.  As  a  descriptor  game  formulation,  the  scalings  are 


confined  to  the  matrix  E .  However,  if  a  state-space  formulation  is  sought,  the  extreme  separation 


of  magnitudes  is  spread  throughout  the  matrices  A  .  B and  B2 ■  This  is  precisely  the  case  which 


degrades  the  numerical  accuracy  of  the  final  results.  Next,  let 


1.2345 

1.98752e  3 

9. be  6 

1.94731 

1 .566e  3 

l.Oe  6 

1  98752e3 

l.Oe  12 

l.Oe  15 

and  S2  = 

1 ,566e  3 

l.Oe  12 

l.Oe  15 

9.0e  6 

l.Oe  15 

1.10102e  38 

l.Oe  6 

l.Oe  15 

1.23976e  35 

The  weighting  matrices  S  [  and  S2  also  contain  elements  with  a  wide  range  of  scalings. 


It  is  customary  to  employ  single  precision  arithmetic  when  demonstrating  numerical 


difficulties.  The  effect  is  to  show  how  serious  errors  can  result  from  even  low-order  problems. 


Assuming  single  precision  (8  significant  digits  of  accuracy),  the  feedbacks  are  computed  using  the 


e  ’ 


1  ■- V  ■’ 
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methods  described  in  the  thesis.  The  result  is 

F ,  =  [  5.06829e  —4  -1.01567e-3  -9.86917e-4] 

F2  =  [-1.73157e -4  1 .22620c  —3  -1.44510e-2]  . 

Now,  multiply  equation  (2.1.1)  by  £“*.  This  yields  three  new  matrices  : 

0.5  -1.00078  -0.99899 

A  =  — 9.88268e  —14  7.77778c -7  -1.0101c -6  , 

9.78385c  —17  -1.95829c -16  -1.94479c -16 

_  '  9.86999c  2  _  1.37000 

Bx  =  1.24223c -6  ,  and  B2=  -1.01037c -9  . 

1.93132c -13  2.68088c -16 

Although  E  -  I.  the  matrices  A.  B  x.  and  B2  are  now  very  poorly  conditioned.  Assuming  the  same 
weighting  matrices,  we  again  compute  the  feedbacks  using  the  best  methods  available.  The  result  is 

Fx  =  [  5.07253e  —4  -1.01652e-3  -9.87737e-4] 

F2  =  [  — 4.78691e  —4  1.84028c-3  -1.3860c -2]  . 

Obviously,  there  is  quite  a  difference  -  especially  in  the  first  element  of  F 2  where  there  is  well  over 

a  factor  of  two  difference.  In  order  to  get  an  accurate  appraisal  of  the  deviations,  consider  the 

descriptor  system  and  run  the  feedback  calculation  using  double  precision  arithmetic  (16  significant 

digits  of  accuracy).  The  result  is 

Fx  =  [5.06600c -4  -1.01521e-3  -9.86498e-4] 

F2  =  [-7.68264c -6  8.94645c -4  -1.47528c -2]  . 

The  reduced  precision  coupled  with  the  scalings  present  in  the  matrix  E  are  responsible  for  the 

decreased  accuracy.  Assuming  this  last  set  of  values  to  be  most  accurate,  we  compute  a  percent 

difference  which  is  summarized  in  Table  2.1. 

The  interesting  item  to  note  for  this  example  is  that  the  state-space  formulation  always  yields 
numbers  which  are  three  times  less  accurate  in  terms  of  percent  difference  from  true  value.  One  is 
led  to  conclude  that  the  descriptor  formulation  is  more  numerically  accurate  (robust)  for  game 
problems  with  nearly  singular  E  matrices.  Clearly,  if  E  is  singular,  the  state-space  formulation  is 
not  even  applicable.  Also,  if  E  is  sparse,  then  multiplying  by  E~ *  is  undesirable  because  sparsity 
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Table  2.1.  Percent  Difference  Calculation 


Component  1 

Component  2 

Component  3 

F  i  Descriptor 

0.045% 

0.004% 

0.042% 

F !  State-Space 

0.13% 

0.13% 

0.125% 

F  2  Descriptor 

2154% 

37% 

2% 

F  2  State-Space 

6131% 

106% 

6% 
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CHAPTER  3 

COMPUTATIONAL  ASPECTS  OF  COUPLED  DISCRETE-TIME  RICCATI  EQUATIONS 

This  chapter  investigates  the  computational  aspects  of  iterating  the  coupled  discrete-time 
Riccati  equations  presented  in  the  last  chapter.  It  was  shown  there  that  coupled  Riccati  equations 
solve  both  single-rate  and  multirate  descriptor  Nash  game  problems.  Though  these  games  have 
been  solved  in  theory,  the  computational  details  of  iterating  coupled  Riccati  equations  are  highly 
nontrivial. 

As  pointed  out  in  Remark  2  following  Fact  2.1,  the  Riccati  equations  are  coupled  through  the 
feedback  matrices.  That  is.  the  ability  to  iterate  the  Riccati  equations  hinges  upon  the  knowledge 
of  both  feedbacks.  But  the  calculation  of  one  DM's  feedback  requires  the  knowledge  of  the  other 
DM's  feedback  (substitute  (2.1.8)  into  (2.1.15)).  This  fact  is  the  source  of  all  numerical  hardships 
encountered  when  attempting  to  iterate  coupled  discrete-time  Riccati  equations.  In  order  to 
produce  any  algorithm  for  iterating  the  Riccati  equations,  the  coupling  must  be  removed. 

Towards  this  end.  Section  3.1  describes  a  feedback  decoupling  procedure.  The  equations  of 
Chapter  2  are  manipulated  in  such  a  way  that  the  coupling  vanishes.  Then  an  algorithm  is  stated 
which  solves  for  the  feedback  matrices  and  iterates  the  Riccati  equations.  Hence  both  Nash  game 
problems  posed  in  the  last  chapter  are  solved  numerically  via  dynamic  programming.  More 
important,  however,  is  the  fact  that  the  task  of  iterating  two  coupled  discrete-time  Riccati 
equations  is  accomplished.  Additionally,  more  notation  is  introduced.  The  section  is  concluded 
with  an  existence  analysis  of  the  algorithm.  Conditions  are  stated  which  insure  that  the  Riccati 
equations  can  be  iterated. 

Section  3.2  presents  the  computationally  superior  coupled  discrete-time  Riccati  iteration 
algorithm.  To  begin  with,  all  relevant  equations  are  assembled  together.  Assuming  the  feedback 
matrices  are  known,  the  problem  of  efficiently  iterating  a  Riccati  equation  is  resolved  by 
considering  numerous  forms  of  the  equation.  Thus  the  final  choice  is  well-justified.  Then  it  is 
observed  that  the  terms  in  the  Riccati  and  feedback  expressions  are  composed  of  positive- 
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(semi)definite.  symmetric  matrices  and  quadratic  forms.  Therefore,  algorithms  are  devised  for 
computing  these  frequently-used  quantities.  Also,  routines  for  multiplying  general  matrices  and 
matrices  with  special  structure  are  given. 

The  task  of  computing  the  feedback  matrices  is  addressed  next.  Unfortunately,  the  (implicit) 
inversion  of  one  general  square  matrix  is  required.  In  fact,  one  of  two  square  matrices  must  be 
inverted.  Hence,  the  computational  burden  is  heavy  at  this  stage  because  the  condition  number  of 
each  matrix  must  be  estimated.  Furthermore,  matrices  with  special  structure  (e.g..  symmetric)  are 
less  prevalent  in  the  feedback  expressions  so  the  computational  intensiveness  problem  is 
exacerbated.  Nevertheless,  it  is  possible  to  rewrite  the  feedback  equations  in  a  form  that  exposes 
additional  structure  thus  alleviating  the  computational  complexity  a  bit.  At  this  juncture,  the 
feedback  algorithm  is  presented.  Finally,  the  coupled  Riccati  algorithm  is  built  from  the  various 
low-level  algorithms  defined  earlier  in  the  discussion. 

Section  3.3  studies  the  existence  of  solutions  to  finite-horizon  problems.  It  is  here  that 
previous  results  (in  Section  3.1)  are  strengthened.  Section  3.4  investigates  the  convergence  behavior 
of  Riccati  iterations  for  infinite-horizon  problems.  The  existence  of  a  region  where  the  coupled 
Riccati  equations  constitute  a  contraction  mapping  is  established.  The  bulk  of  the  section  is  devoted 
to  proving  the  contraction.  The  contraction  mapping  argument  guarantees  existence  of  and 
convergence  to  a  unique  fixed  point  for  an  infinite  number  of  Riccati  iterations. 

3. 1  A  Decoupling  Procedure 

This  section  considers  some  preliminary  computational  aspects  of  iterating  the  coupled 
discrete-time  Riccati  equations  obtained  in  the  previous  chapter.  First,  the  Riccati  equations  are 
decoupled.  Then  an  LQ  Nash  game  algorithm  is  presented.  Finally,  various  necessary  and  sufficient 
conditions  which  govern  existence  and  uniqueness  of  solutions  to  finite-horizon  LQ  Nash  game 
problems  (and  hence  coupled  Riccati  iterations)  are  stated. 
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Beginning  with  (2.1.6)  use  (2.1.8).  (2.1.10)  and  (2.1.13)  to  obtain  : 
u  j*(k)  =  -F1(k)A:,(k)  =  -(r^k))"1  [5[A:,(k+l)  {A  -B2F2{ k))]  x*(k) 


(3.1.1) 


u2*(k)  =  -F,(k)x‘(k)  =  -(r2(k))"1  [B$K2(k+l)  ( A  -B^J k))]  x'(k)  .  (3.1.2) 

Notice  that  (3.1.1  )-(3. 1.2)  are  coupled  equations  where  the  coupling  occurs  through  the 
feedback  matrices  Fj(k)  and  F2(k)  .  In  order  to  present  an  algorithm  for  solving  the  Riccati 
iteration  problem,  this  coupling  must  be  addressed.  Towards  this  end.  we  will  use  a  cross¬ 
substitution  procedure  to  derive  the  functional  relationships  : 

F,(k)  =  g,U  jfk+l).  tf2(k+l))  (3.1.3a) 

F2(k)  =  g2U,(k  +  l).  A'2(k+1))  .  (3.1.3b) 

To  simplify  notation,  we  introduce 

V,(k)  i  (r.(k))-1  B\  A\(k+1).  i=l .2.  (3.1.4) 

Therefore,  from  (3.1.1)  and  (3.1.2)  we  conclude  that  : 

F,(k)  =  (T^k))"1  B\K  j(k+l )  (A  -B2F2( k))  (3.1.5a) 

=  *,(kM  -  *,(k)52F2(k)  (3.1.5b) 

=  V,(k).4  -  ♦^k^z^T^k))"1  5fA-2(k+l)  (.4  -fl,F,(k))  (3.1.5c) 

=  *,(k)  I-52^2(k)  .4  +  ^1(k)fl2^2(k)J5,F1(k)  .  (3.1.5d) 

Hence. 


1  -  'k1(k)/?2'J'2(k)fl, 

i  C  k )  =  ¥,(k) 

I  -  B  2¥2(  k ) 

Similarly,  for  the  feedback. 

/•'2(k).  we  get  : 

I  -  V2(k)fl,V,(k)fl2 

F2(  k)  =  V2(k) 

1  -  A,Y,(k) 

Define. 


=  ,(k)  “  1  -  ¥1(k)fl2*2(k)/?I 


E2(k)  fj  1  -  ♦2(k)«1'I'1(k)«2 


(3.1.4) 

K/ 

K 

►  * 

(3.1.5a) 

P 

(3.1.5b) 

(3.1.5c) 

*V 

(3.1.5d) 

y 

(3.1.6a) 

% 

(3.1.6b) 

>. 

V 

(3.1.7a) 

■ 

(3.1.7b) 

.* 

nwww 


Then  the  functions  g  i  and  g2  described  in  (3.1.3)  are  given  by 


I! 

FJ  k)  = 

(si(t))"1 

a 

F2{  k)  = 

(s.ik)!-1 

\  / 

> 

assuming  that 

l^(k))-1 

V 

1 

_ 1-1 

(3.1.8a) 

(3.1.8b) 


bounded  entries.  Existence  of  |£i(k)J  i=1.2  is  a  more  delicate  issue  that  will  be  dealt  with 

shortly. 

We  may  now  present  an  algorithm  for  solving  the  LQ  descriptor  Nash  game  problem. 
Dynamic  programming  [5]  dictates  that  the  feedbacks  must  be  solved  for  in  reverse  time.  The 
requirement  of  iterating  two  coupled  discrete-time  Riccati  equations  is  automatically  satisfied. 

IjQ  Descriptor  Nash  Game  Algorithm 


Step  1  :  Initializations 
Set  COUNT  -  T,  . 

Set  Er  K .(COUNT  +  1)  E  =  Cf  S .{COUNT  +  1 )  C,.  i=l  .2. 


(Terminal  Constraint) 


Step  2  :  Beginning  of  Main  Loop 

Compute  T .{COUNT)  .  i= 1,2  using  (2.1.13). 

Step  3  :  Core  Calculations 

Compute  'if .{COU NT  ),  i=1.2  using  (3.1.4) 

Compute  E. {COUNT  ),  i=1.2  using  (3.1.7). 

If  \~:{COUNT )  exists,  then  compute  it.  otherwise  go  to  Step  5. 

If  ^E itCOU NT  )  exists,  then  compute  it.  otherwise  go  to  Step  5. 

Compute  F .{COUNT  ).  i=l  .2  using  (3.1.8). 

Compute  A  .{COUNT  ).  i=1.2  using  (2.1.8). 

Step  4  :  Update  Data 

Store  F .(COUNT  ).  i=1.2  for  this  value  of  COUNT  . 

Iterate  the  coupled  discrete-time  Riccati  equations  according  to  (2.1.11)  with  i  =  1 .2  thus  obtaining 
K. {COUNT  )  from  A  .(COUNT  +1 ).  i=  1 .2 
Set  COUNT  -  COUNT  -  1 

If  COUNT  ^  0.  then  go  to  Step  2.  otherwise  STOP  the  algorithm. 


**  V  *.  V 
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Step  5  :  Error  Condition 

Print  an  Error  Message  that  indicates  that  one  or  more  of  the  feedbacks  are  infinite  at  this  stage  of 
the  game.  Then  STOP  the  algorithm. 

Remark  :  Multirate  LQ  Descriptor  Nash  Game  Algorithm 

The  corresponding  multirate  algorithm  may  be  obtained  by  simply  replacing  step  2  with 
Compute  Ir^GOWDr1.  i— 1.2  using  (2.1.13).  If  L*Z+,  then  Set  \t2(COUNT  )  I  *  =  0. 


Moreover,  other  simplifications  are  possible  for  a  multirate  game.  If  L  (f  Z+.  then 
|r2(C06Wnj  =  0  which  implies  that  V2{COUNT )  =  0.  Furthermore.  (3.1.7)  implies  that 

S i(COUNT )  =  £2(COUNT )  =  I  which,  in  turn,  implies  that  F^COUNT)  =  Vl(COUNT)A 
and  F2(COUNT  )  =  0.  Hence,  the  Riccati  equations  that  have  to  be  iterated  reduce  to 

Er  Kt(k)E  =  AT  Jfi(k+1)- A1(k+l)31(ri(k))-1flfA1(k+l)  A  +  G^k) 

ET  A'2(k)  E  =  A 2  (k)A’2(k+l)A 2(k)  +  fi2(k)  . 

Many  facts  become  apparent  from  the  descriptor  Nash  game  algorithm.  They  are  summarized 
by  the  following  : 

Proposition  3.1  :  Existence  of  Feedbacks 

Aj(k)  exists  if  and  only  if  I  fSjfk)  1  I  <  oo  Likewise.  A2(k)  exists  if  and  only  if 


II  |H2(k)  fl  <  oo  .  Equivalently,  there  exists  a  unique  solution  (y*  .  y2k  )  ,  k  €  H  to  the  kf/l  stage 

of  the  LQ  descriptor  Nash  game  described  in  Chapter  2  if  and  only  if  crfc^k))  &  0  and 
(r(s2(k))  ^  0  for  that  k.  Therefore,  there  exists  a  unique  solution  policy  (yj,  y2)  to  the  entire 
Nash  game  if  and  only  if  0  <E  {  a(Hs(k)  I  k€H,  i=  1 .2 } . 

Proof  :  Consider  only  single-rate  games  and  let  k =T f  .  Then.  A,(k+l)  .  i=1.2  are  well-defined 
by  (2.1.14).  Assumption  2.1  guarantees  that  A',(k+1)  is  unique.  Since.  0  <  A  j  <  oo  and 
II  B\  II  <  oo,  T[(k)  exists.  Similarly.  0  <  R2  <  °°  and  II  B2  II  <  oo  implies  that  T2(k)  exists. 

Thus.  [r,(k)|  ,  i=  1 .2  exists  if  and  only  if  A,(k+1)  exists,  respectively. 
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Now,  given  that  Jr^k)  J  ,  i— 1.2  exists,  then  by  (3.1.4)  so  does  ^(k) .  Furthermore,  from 

(3.1.8),  it  is  clear  that  F’j(k).  i=l,2  exists  if  and  only  if  the  corresponding  |Sj(k)j  exists.  This 

condition  is  the  same  as  requiring  <r(Ej(k))  0  ,  i«1.2.  From  (2.1.11a),  it  is  obvious  that 

boundedness  of  A"i(k)  hinges  upon  boundedness  of  Aj(k)  which,  from  (2.1.8),  occurs  if  and  only  if 

I  Fj(k)  I  <  oo.  i=l,2  .  Therefore.  X'j(k)  exists  if  and  only  if  a  [(S^k))-*  |  = - - -  <  oo  .  That 

1  '  a(S,(k)) 

is.  q[(sj(k))  ^  0.  i=1.2.  Uniqueness  comes  from  the  fact  that  (2.1.11a)  defines  a  single  matrix, 
ET  if  i(k)  E  given  all  quantities  on  the  right-hand  side  of  (2.1.11a)  and  Assumption  2.1.  Finally, 
the  proof  is  completed  by  inductively  applying  the  above  argument  to  the  next  (i.e..  k— 1”  )  stage 
of  the  game  as  long  as  {if^k).  i=l  .2}  exists.  If  not.  then  the  given  LQ  descriptor  Nash  game  is  ill- 
posed  because  either  u  j  (k)=oo  or  u\  (k)=oo  or  both  for  some  k  €  H. 

□ 


Remark  :  Existence  of  Multirate  Feedbacks 

If  L  €  Z+,  then  the  result  of  Proposition  3.1  applies.  However,  whenever  L  t  Z+,  the  remark 
following  the  LQ  descriptor  Nash  game  algorithm  indicates  that  both  feedbacks  always  exist. 
Notice  that  if  N  =2,  then  the  conditions  L  €  Z+  and  L  t  Z+  occur  with  equal  frequency.  But  if 
N  >2,  then  L  i  Z+  occurs  more  often.  Hence,  as  N  increases,  the  frequency  with  which  the 
feedbacks  are  guaranteed  to  exist  increases. 

As  an  immediate  consequence  of  Proposition  3.1,  we  have  : 

Corollary  3.1  :  A  sufficient  condition  for  existence  of  a  unique  solution  to  the  finite-horizon  LQ 
.Nash  game  problem  is  that  (ff’PjCk)/?  2¥2(k)Z?  j)  <  1  and  <  1  for  all  k  €  H. 


i 

•» . 
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Proof  It  is  well  known  [18.47]  that  the  singular  values  of  a  matrix  are  intimately  related  to  the 
problem  of  rank  degeneracy.  In  particular,  given  a  square  matrix  X  with  cr(X  )  ^  try  >0  and 

another  square  matrix  Y  of  compatible  dimension  with  cf(l' )  ^  ^0.  then  if  cf5  <  <r.v  the 


matrix  (X  +Y )  has  full  rank.  Consider  (3.1.7a)  and  take  X  =  I  and  Y  —  ¥i(k).ff2*2(k)Zfi  • 
Obviously,  <f(I)  =  gi(I)  =  1.  Hence  as  long  as  <7(Y )  =  cF^^k)/?^^)/? 1)  <  1.  the  matrix  Ej(k) 
can  never  be  singular.  Thus.  [Si(k)  I  always  exists.  Similarly,  consider  (3.1.7b)  and  take 


Y  =  ^(k)/?  i'ki(k)52  leaving  X  as  before.  The  same  singular  value  argument  can  be  applied  to 
this  case,  which  then  proves  the  Corollary. 


Corollary  3.2  :  Another  sufficient  condition  for  existence  of  a  unique  solution  to  the  finite-horizon 
LQ  Nash  game  problem  is  that  tr^ifk^^^k)/?  i)  >  1  and  gtft^k^i^fk)^)  >  *  ^or  ^  ® 

H. 

Proof  :  Using  the  same  singular  value  argument  as  discussed  in  Corollary  3.1.  consider  (3.1.7a) 
and  take  Y  =  I  and  X  =  Yi(k)£2¥2(k)£i  -  As  long  as  <Zl(X  )  >  <F(K )  =  1.  the  matrix  Ej(k)  can 
never  be  singular.  Similarly,  consider  (3.1.7b)  and  take  X  =  ^2(k)5 1^,(k)^2  leaving  Y  as  before. 
The  same  singular  value  argument  can  be  applied  to  this  case,  which  then  proves  the  Corollary. 


We  can  carry  these  last  two  corollaries  one  step  further  by  restricting  attention  to  those 
games  in  which  both  players  have  a  single  input  to  the  system.  In  this  case,  define  the  scalars  : 

w,(k)  *  'MkjB, 
u3(k)  ^  Y,(k)B,  . 

Then.  Ei(k)  =  1  —  <t>)(k)oj2(k)  =  1  —  (ujCkltu^k)  =  E2(k)  .  Hence, 

Corollary  3.3  :  A  unique  solution  to  the  finite  horizon  I.Q  Nash  game  problem  with  single  inputs 
for  each  player  exists  if  and  only  if  a>j( k )  ^  (aj2(k))-*  for  all  k  €  H. 

Proof  Obviously.  E((k)  =  E2(k)  are  scalar  quantities  for  the  situation  when  each  player  has  a 
single  input.  Therefore,  singularity  of  both  occurs  if  and  only  if 
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<«»1(k)ft>2(k)  =  fc^Ck^^k)  =  1 
This  condition  is  equivalent  to 

G>i(k)  =  (^(k))-1  or  ^(k)  =  (^(k))-1  . 

□ 

In  summary,  these  results  provide  on-line  checkable  conditions  so  that  the  descriptor-game 
algorithm  can  monitor  itself  and  determine  if  a  Nash  solution  strategy  exists.  In  the  next  section, 
algorithms  are  developed  for  iterating  coupled  discrete-time  Riccati  equations  in  a  computationally 
efficient  manner. 

3.2  The  Coupled  Discrete— Time  Riccati  Algorithm 

In  Section  3.1  an  algorithm  is  presented  for  solving  LQ  descriptor  Nash  games.  The  process 
involves  iterating  coupled  discrete-time  Riccati  equations.  This  section  describes  a  more  efficient 
and  numerically  robust  procedure  for  performing  one  iteration  of  two  coupled,  discrete-time 
Riccati  equations.  As  a  byproduct  of  the  analysis,  we  obtain  an  equally  efficient  procedure  for 
iterating  a  single  discrete-time  Riccati  equation  which  arises  in  the  study  of  the  optimal  LQ 
regulator  problem. 

3.2.1  Algorithm  Implementation 

This  subsection  reviews  the  pertinent  equations  involved  in  iterating  two  coupled  discrete¬ 
time  Riccati  equations.  A  complete  discussion  of  the  equations  may  be  found  in  Section  3.1. 
Because  positive~(semi)definite  symmetric  matrices  are  frequently  encountered  in  this  coupled 
Riccati  problem,  the  Cholesky  factorization  [17.24]  will  be  used  to  expedite  the  calculations.  Recall 
that  any  positive-(semi)definite  symmetric  matrix,  A  .  may  be  factored  as  A  —  XAXA  where  XA  is 
upper  triangular.  It  is  preferable  to  work  with  X A  rather  than  directly  with  A  because  of  the 
triangular  structure.  For  example,  it  is  often  necessary  to  calculate  C  =  A~^  B  where  B  is  some 
other  compatibly-dimensioned  matrix.  This  can  be  accomplished  by  solving  the  symmetric  linear 
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system  A  C  =  B  for  C .  However,  if  the  Cholesky  factor  of  A  is  available  then  the  linear  system 
is  equivalent  to  solving  two  triangular  systems  XA  D  =  B  and  XA  C  =  D .  This  technique  has 
two  major  advantages  over  working  directly  with  the  original  system  A  C  —  B .  First,  triangular 
systems  are  extremely  easy  to  solve.  Second,  the  inverse  of  A  is  never  explicitly  computed.  As  a 
general  rule  of  thumb,  the  explicit  computation  of  an  inverse  should  be  avoided  at  almost  any  cost. 
For  the  case  where  A  is  dense  and  has  no  special  structure,  the  Cholesky  factorization  is  replaced 
by  the  LU  factorization.  In  this  factorization.  U  is  upper  triangular  and  L  is  the  product  of 
elementary  lower  triangular  and  permutation  matrices. 

Remarks  :  Sof tware  Engineering 

The  solution  of  this  coupled  Riccati  problem  requires  software  designed  to  handle  linear  algebraic 
equations.  Frequently.  positive-(semi)definite  symmetric  matrices  appear  which  require  special 
handling.  The  UNPACK  [24]  FORTRAN  library  is  the  best  known  and  the  most  widely 
recommended  [21]  software  package  used  for  these  situations.  For  all  the  algorithms  which  follow. 
Single  (S)  precision  arithmetic  is  assumed.  Therefore,  only  subroutines  beginning  with  'S'  are 
referenced  from  the  UNPACK  library.  If  Double  (D)  precision  is  used,  then  the  first  letter  of  the 
corresponding  LINPACK  routine  is  D'.  However.  L-A-S  is  a  double  precision  package.  Since  the 
coupled  Riccati  software  will  be  integrated  into  L-A-S.  the  software  listings  in  Appendix  A  refer  to 
the  double  precision  versions  of  LINPACK.  Second,  all  of  the  Riccati  software  discussed  here  is 
coded  in  FORTRAN.  This  decision  is  primarily  motivated  by  the  fact  that  both  L-A-S  and 
LINPACK  are  FORTRAN-based  packages.  In  addition.  FORTRAN  is  generally  chosen  for  coding 
numerical  software. 

Several  of  the  algorithms  require  at  least  one  multiplication  of  two  matrices  with  no  special 
structure  (like  upper/lower  triangular).  For  this  reason,  we  define  algorithm  .MLTPLY,  which 
multiplies  two  arbitrary,  but  compatible  matrices.  In  particular  given  the  matrices  A  and  B  , 
MLTPLY  is  useful  for  computing  the  forms  A  B  and  A  '  B  .  Notice  that  the  case  where  both  A 
and  B  are  symmetric  offers  no  advantages  since  the  resulting  matrix  is.  in  general,  not  symmetric. 


>  ''oHI 


Hence,  algorithm  MLTPLY  is  applicable.  Since  MLTPLY  is  too  simple  to  be  given  as  a  sequence  of 
steps,  a  copy  of  the  FORTRAN  source  code  is  included  in  Appendix  A. 

Now  the  equations  needed  to  iterate  coupled  Riccatis  are  reviewed.  The  Riccati  equation  for 
Decision  Maker  i  (DMi).  i  €  N  ^  {  1.  2}  at  stage  k  is  given  by 

ET  AT t(k)  E 


=  Af  (k)A',(k+l)A1(k)  -  Ar(k)A’l(k+l)5,(ri(k))-1  BfA-,(k+l)A  >(k)  +  Q,(k)  (3.2.1) 


where 


r,(k)  2  *,(k)  +  BfK,(k+l)B,  . 
A  ,(k)  J  A  -  B2F2(k )  .and 
A  2(k)  *  A  -fl,F,(k). 


The  feedback  employed  by  each  DM  is 

=  (/?  ,(k)  +  B\K  y{k  +  \)B  (/?(A,(k+l)A  ,(k)) 

=  (T^k))-1  [B(Ar,(k  +  l)  (A  -B2F2( k))] 

=  |H,(k)|_1  ¥,(k)  |  I  -  i?2*2(k)|  A 

and 

F2(k)  =  (/?2(k)  +  ^/^(k+l)^)"1  (fi(A2(k+l)A2(k)) 
=  (T^k))"1  (5CA-2(k+l)  (A  -#,F,(k))] 


=  =2( k )  ¥2(k)  I-fli^i(k)  .4 


where 


♦  ,(k)  “  (r.(k))-1  B{  A,(k+1).  i=1.2 


H,(k)  *  I  -  V,(k)J92¥2(k)2?,  .  and 
=2(k)  ^  1  -  ♦2(k)/}1'k,(k)«2  . 


(3.2.2) 


(3.2.3a) 


(3.2.3b) 


(3.2.4a) 

(3.2.4b) 


(3.2.4c) 


(3.2.5a) 

(3.2.5b) 


(3.2.5c) 


(3.2.6) 


(3.2.7a) 


(3.2.7b) 
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3.2.2  Iterating  a  Riccati  Equation 


This  subsection  determines  an  algorithm  for  iterating  a  Riccati  equation  assuming  the 


feedbacks  are  known.  The  goal  here  is  to  minimize  the  number  of  inverses  appearing  in  equation 


(3.2.1).  From  a  computational  viewpoint,  the  right-hand  side  of  the  Riccati  equation  (3.2.1)  is 


most  efficiently  calculated  by  using  a  slightly  modified  form  of  (3.2.1).  To  see  this,  start  with 


(3.2.1)  and  use  the  relations  (3.2.4)  and  (3.2.5)  to  write 


ET  K,(k)E 


=  Af  AT,(k+l)  -  Af  XiCk+DBiCr,)-*  Bf  K  ,(k +  1M,  +  Q, 


(3.2.8a) 


=  A  f  Ktk+DA,  -  A[  Kfk+DBiCr,)-*  f,  ff 1  Bf  A,(k+Di4,  +  Q ,  (3.2.8b) 


=  Af  tf.U+lM,  -  (Af  AT,(k+l)  B,  Bf  K  j(k+l)  A ,)  +  <2, 


(3.2.8c) 


=  Af  A.U+lM,  “  Ff  r ,F,  +  Q , 


(3.2.8d) 


where  the  dependence  on  k  has  been  dropped  for  all  variables  except  K 


Remark  :  It  is  a  simple  exercise  to  verify  that  (3.2.8d)  can  be  written  as: 


ET  A'.(k)  E  =  A{l  A,(k+l).4ri  +  Ff  /?,A,  +  <2, 


(3.2.8e) 


where  ACL  _  A  —  BlFl~B2F2-  This  form  is  as  compact  as  (3.2.8d).  Indeed,  this 


computationally  attractive  form  may  be  found  in  [7,p.253]  where  the  setting  is  largely  theoretical. 


However,  it  will  become  evident  that  it  is  less  efficient  to  use  (3.2.8e)  to  iterate  the  Riccati 


equations,  because  the  Cholesky  factor  of  /?,  is  never  computed  whereas  the  Cholesky  factor  of  T, 


will  be  available  from  another  calculation. 


One  additional  simplification  is  possible.  In  view  of  (3.2.6).  equation  (3.2.8c)  can  be  rewritten 


A  A  ,( k )  /:  =  ,1;  A',(k+1)  -  f,  .4,  +  Q, 


( 3.2.M  ) 


which  is  extremely  compact. 


It  is  true  that  each  Riccati  equation  can  be  computed  directly  as  a  (unction  ol  A'j(k)  anti 


Ai(k).  but  the  resulting  equation  is  very  complex  ami  quite  ol  ten  the  feedbacks  are  needed 


,  *  .  4 
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Furthermore,  once  the  feedbacks  are  calculated,  the  Riccati  iterations  are  most  efficiently  computed 
using  (3.2.8f)  -  an  equation  with  no  inverses. 

Before  presenting  the  algorithm  for  iterating  a  Riccati  equation,  several  small  algorithms  must 
be  defined.  This  is  done  in  the  next  few  subsections.  Algorithm  QDFORM  computes  a  quadratic 
form  or  variations  of  a  quadratic  form.  Algorithm  PSICOM  computes  the  variable  ’Pj(k)  as  defined 
in  (3.2.6).  Algorithm  XTRACT  computes  K ,(k)  given  the  right-hand  side  of  (3.2.8f).  Algorithm 
DISR1C  combines  these  algorithms  to  perform  one  iteration  of  a  Riccati  equation. 


3.2.2. 1  Algorithm  QDFORM 

It  is  clear  that  the  quadratic  form  BT  A  B  where  B  is  arbitrary  and  A  is  a  positive- 
(semi)definite.  symmetric  matrix  appears  frequently.  Anticipating  the  need  to  compute  several 
quadratic  forms,  let  us  introduce  the  following  algorithm. 


Algorithm  QDFORM 


Step  1  :  Factor  A  _  XA  XA  where  XA  is  upper  triangular.  XA  is  called  the  Cholesky  factor  of  A 

and  is  obtained  via  the  LINPACK  subroutine  SPOFA. 

Step  2  :  Compute  Y  =  XA  B  taking  advantage  of  the  fact  that  XA  is  upper  triangular. 

Step  3  :  Compute  D  -  Yr  Y  taking  advantage  of  the  fact  that  D  is  symmetric. 

Remarks  :  As  coded  in  Appendix  A.  Step  1  can  be  bypassed  if  the  Cholesky  factor  of  A  is  already 
available.  Steps  2  and  3  require  one  optimized  matrix  multiply  DO-Loop  in  FORTRAN.  Also,  if 
the  quantity  C  +  BT  A  B .  where  C  is  symmetric,  is  desired,  then  algorithm  QDFORM  is 
applicable  if  Step  3  is  modified  to  include  the  array  C  as  the  initial  condition  in  the  multiply  DO- 
Loop.  The  same  can  be  said  for  the  form  C  —  BT  A  B  . 


42 


3.2.2.2  Algorithm  PSICOM 

The  computation  of  'Pj(k),  i=1.2  establishes  the  foundation  upon  which  all  subsequent 
calculations  are  built.  Observe  that  (3.2.4),  (3.2.5).  (3.2.7)  and  (3.2.8)  depend  on  ¥j(k).  It  is 
apparent  that  a  thorough  investigation  of  the  computation  of  ¥j(k)  is  needed.  From  the  definition 
(3.2.6),  it  would  seem  that  the  explicit  computation  of  (r,(k))-^  is  necessary.  However,  this  is  not 
the  case.  Since  rs(k)  is  a  positive-definite,  symmetric  matrix,  the  LINPACK  [24,  Chapter  3] 
software  for  general  positive-definite  matrices  is  used  to  circumvent  this  problem. 

Although  ¥j(k)  can  be  found  in  4  easy  steps,  the  result  of  each  step  of  the  calculation  is 
needed  in  later  computations.  For  example,  the  quantity  W  defined  in  Step  1  of  algorithm  PSICOM 
is  used  later,  so  let 

W,( k)  *  B\  K  1(k+l).  and  (3.2.9a) 


W,(k)  *  B\  K 2(k+l )  .  (3.2.9b) 

Hence  it  is  desirable  to  create  an  algorithm  that  computes  and  returns  ¥j(k)  as  well  as  all  the 
intermediate  quantities.  Since  4^,(k)  is  a  function  of  Tj,  B,  and  K t.  the  subscript  i  can  be  dropped. 
Also,  the  dependency  on  k  can  be  suppressed.  The  algorithm  for  computing  'It  is  now  stated. 

Algorithm  PSICOM 

Step  1 :  Compute  W  ^  Br  K  using  algorithm  MLTPLY. 

Step  2  :  Compute  T  using  algorithm  QDFORM. 

Step  3  :  Factor  T  =  X[  Xr  where  Xr  is  upper  triangular.  Xv  is  called  the  Cholesky  factor  ot  T 
and  is  obtained  using  the  LINPACK  subroutine  SPOFA. 

Step  4  :  Form  'P  —  I  ^  W  without  explicitly  computing  using  the  LINPACK  subroutine 
SPOSL. 


*  m 


& 


V 
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Remarks  :  Step  4  solves  the  linear  equation  T  'P  =  W  for  'P  one  column  at  a  time  using  the 
Cholesky  factor  Xr.  Thus,  it  is  not  necessary  to  explicitly  form  the  inverse  of  f  Algorithm 
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PSICOM  takes  B .  K .  and  £  as  input  and  subsequently  produces  W .  T,  X r.  and  'If  as  output. 


Hence,  if  'If  is  needed  for  a  later  calculation  and  algorithm  PSICOM  is  invoked,  then  f  need  not  be 


computed  using  algorithm  QDFORM  because  PSICOM  will  return  it  as  a  byproduct  of  the 


calculation  of  ¥.  The  Riccati  algorithms  will  exploit  this  savings  in  computation. 


3.2.2.3  Algorithm  XTRACT 


Given  the  right-hand  side  of  (3.2. 8f),  the  quantity  ATj(k)  is  immediately  known  only  if  E  —  I. 


If  E  differs  from  the  identity  matrix,  then  additional  steps  must  be  taken  to  extract  ATjCk)  from 


the  quadratic  form  ET  K  fX)  E .  Algorithm  XTRACT  performs  exactly  that  function.  Specifically. 


suppose  the  right-hand  side  of  (3.2.8f)  evaluates  to  a  positive-definite,  symmetric  matrix  called  Z  . 


Then  Z  has  a  Cholesky  factorization.  Z  _  Xz  Xz  .  Similarly,  the  kernel  of  the  left-hand  side  of 


(3.2. 8f)  has  a  Cholesky  factorization  K  ^  X%  XE  .  Equating  both  quantities  yields 


ET  X[  XK  E  =  X{  Xz  . 
Algorithm  XTRACT  is  based  on  this  observation. 


Since  E  is  assumed  to  be  dense  and  with  no  special  structure,  the  LINPACK  software  for 


general  matrices  is  employed.  Assumption  2.1  guarantees  that  the  inverse  of  E  exists.  However,  as 


with  algorithm  PSICOM.  the  inverse  is  never  explicitly  computed.  Instead  the  LU  factorization  of 


E  is  used.  Now.  the  algorithm  for  extracting  A",(k)  is  stated. 


Algorithm  XTRACT 


Step  1  :  Factor  Z  -  Xz  X7  where  Xz  is  upper  triangular.  Xz  is  called  the  Cholesky  factor  of  Z 


and  is  obtained  using  the  LINPACK  subroutine  SPOFA. 


Step  2:  Factor  E  -  LE  UE  where  UE  is  upper  triangular  and  LE  is  the  product  of  elementary 


lower  triangular  and  permutation  matrices.  LE  UE  is  called  the  LU  factorization  of  E 


and  is  obtained  using  the  LINPACK  subroutine  SGEFA. 
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Step  3:  Solve  the  linear  system  Er  Xg  =  X}  for  XK  without  explicitly  computing  E~ *  using 
the  LINPACK  subroutine  SGESL. 

Step  4 :  Compute  K  =  Xf  Xg  taking  advantage  of  the  fact  that  D  is  symmetric. 

Remarks  :  Step  3  solves  the  linear  equation  ET  Xg  =  X[  for  XK  one  column  at  a  time  using  the 
UJ  factorization  LE  UE .  Thus,  it  is  not  necessary  to  explicitly  form  the  inverse  of  E.  Step  4 
requires  one  optimized  matrix  multiply  DO-Loop  in  FORTRAN.  Algorithm  XTRACT  is  only 
invoked  if  E  &  I.  However,  if  E  should  possess  additional  structure  (e.g..  diagonal),  then  steps  2 
and  3  should  be  replaced  by  another  algorithm  which  exploits  that  structure. 

3.2.2.4  Algorithm  DISR1C 

This  algorithm  combines  the  previous  algorithms  to  perform  one  iteration  of  a  Riccati 
equation.  It  is  assumed  that  the  feedbacks  (  and  therefore  A ,.  i€N  )  are  known.  First  the  right- 
hand  side  of  (3.2.8f)  is  computed.  Then  A'/k)  is  extracted.  The  algorithm  for  iterating  a  Riccati 
equation  is  stated  below. 

Algorithm  DISR1C 

Step  1 :  Compute  ¥  using  algorithm  PSICOM.  T  and  Xv  will  also  be  computed  and  provided  upon 
return. 

Step  2 :  Form  V  =  K  —  Vr  T  V  using  algorithm  QDFORM  taking  advantage  of  the  fact  that 
Xr  is  already  available. 

Step  3  :  Form  Z  =  Q  +  Ar  V  A  using  algorithm  QDFORM. 

Step  4 :  Compute  A'  C k )  from  the  right-hand  side  of  (3.2.8f)  using  algorithm  XTRACT.  if 
necessary. 


Remarks  :  Notice  that  every  step  of  algorithm  DISRIC  is  a  call  to  one  of  the  previously-defined 
algorithms.  This  indicates  that  the  low-level  routines  are  very  modular.  This  is  one  ol  the 
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trademarks  of  a  good  structured-programming  approach. 


3.2.3  Calculating  the  Feedbacks 

Having  devised  a  method  for  iterating  a  Riccati  equation,  we  turn  our  attention  to  calculating 
the  feedbacks.  Using  (3.2.6)  we  can  rewrite  (3.2.4a)  and  (3.2.5a)  as 


f\(k)  =  ^(kMjOt).  i  €  N.  (3.2.10) 

But.  this  form  assumes  knowledge  of  the  other  DM's  feedback  matrix.  Initially  neither  F  j  nor  F 2 

is  known.  Hence  either  (3.2.4c)  or  (3.2.5c)  must  be  used  to  compute  one  of  these  matrices  based 

upon  some  criterion.  Each  equation  is  undesirable  because  it  contains  an  inverse.  Thus,  it  is 

ultimately  necessary  to  compute  the  inverse  of  one  general  square  matrix.  Furthermore,  the 

decoupling  procedure  used  to  produce  (3.2.4c)  and  (3.2.5c)  clouds  the  presence  of  symmetric  and 

positive-definite  matrices.  Therefore,  it  is  worthwhile  to  study  these  equations  with  the  intention 

of  exposing  any  additional  symmetry  and/or  positive  definiteness.  Towards  this  end  take  (3.2.4c) 


and  premultiply  by  H^k).  The  result  is 

s,(k)/’,(k)  =  ¥,(k)  I-fi2¥2(k)  A  . 

Substitute  (3.2.6)  and  (3.2.7a)  into  (3.2.11a)  to  obtain 

1 1  -  (r,(k))“1fi[/i:1(k+l)52(r2(k))-1j5jA'2(k+l)51 )  F,(k) 

=  (T1(k)r15[/irl(k+l)/l  -(r1(k))-15(A'1(k+l)52(r2(k))-1fi5A:2(k+lM 

Premultiplying  by  Tj(k)  yields 

r,(k)-5[7ir1(k+i)52(r2(k))-15?/i:2(k+i)51)  F,(k) 

=  B[K  [(k+l  X4  -5fA'1(k+l)fi2(r2(k))-15^A'2(k+lM  . 

Define  the  positive-(semi)definite,  symmetric  matrices 

X,(k)  t  (r.(k))-1  B\  .  i  €  N  . 

Then  in  view  of  (3.2.9)  and  (3.2.12).  equation  (3.2.1  lc)  reduces  to 


(3.2.11a) 


(3.2.11b) 


(3.2.11c) 


(3.2.12) 


T,(k)  -  W1(k)X2(k)/f2(k+l)51  /■-,( k)  =  W,(k)  1  -X>(k)  A'2(k+1)  .4  .  (3.2.13a) 


A  similar  derivation  applied  to  (3.2.5c)  gives 

|r2(k)  -  W2(k)Xi(k)A'1(k+l)52|  f2(k)  =  W2(k)  1 1  —  Xi(k)  A’lCk+l)  J  A  .  (3.2.13b) 

Define 

Yj(k)  *  r,(k)  -  W x(k)  X2(k)  X2(k+l)  B  ,.  and  (3.2.14a) 

Y2(k)  *  r2(k)  -  VT2(k)Xi(k)X,(k+l)52  .  (3.2.14b) 

Then,  it  is  wise  to  choose  Yj(k)  such  that  it  has  smallest  condition  number  over  all  other  Y/s. 
i  €  N.  Specifically,  suppose  that  DM  i  has  Y-(k)  with  smallest  condition  number.  Obviously  we 
want  to  compute  j  Yj  J  first. 

It  is  apparent  that  the  matrices  Xi(k),  i  €  N  are  needed  for  the  feedback  calculation.  Since 

algorithm  QDFORM  is  not  suited  to  handle  an  inverse  as  the  kernel  of  the  quadratic  form, 

algorithm  CHICOM  is  presented. 

Algorithm  CHICOM 

I  |-7' 

Step  1 :  Compute!  =  Xr  Br  using  UNPACK  subroutine  STRSL.  Xr  is  the  Cholesky  factor 
of  T  obtained  from  algorithm  PSICOM. 

Step  2  :  Compute  X  =  YT  Y  taking  advantage  of  the  fact  that  X  is  symmetric. 

The  following  algorithm  determines  i  and  also  computes  Y,(k),  for  all  i  €  N.  Since  the  final  result 
of  this  algorithm  is  the  computation  of  i.  it  is  most  advantageous  to  code  this  as  an  INTEGER 
FUNCTION  in  FORTRAN  which  returns  an  integer  equal  to  i.  Note  that  for  this  algorithm  the 
dependence  on  i€  N  cannot  be  suppressed. 

Algorithm  EYEHAT 

Step  1  :  Compute  Xi(k).  i  6  N  using  algorithm  CHICOM. 

Step  2:  Compute  7*,C k)  =  VV  ,(k )  X2U )  and  7N(k)  =  W,(k)Xi(k)  using  algorithm  MI.TPLY. 
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Step  3  %  Compute  V  j(k)  =  T  ,(k)  A2(k+l)  and  V2(k)  =  T 2(k)  A  ,(k+l)  using  algorithm  MLTPLY. 
Step  4:  Compute  Yi(k)  =  T^k)  —  Vj(k)  5j  and  Y2(k)  =  T2(k)  —  V2(k)  B2  in  one  matrix 
multiply  DO-Loop. 

Step  5  s  Factor  Yj(k)  =  L  Y  f/Y,  and  Y2(k)  =  L  Yj  i/Yj  where  t/Y.  is  an  upper  triangular  matrix  and 
Z.Y.  is  the  product  of  elementary  lower  triangular  and  permutation  matrices,  i  €  N. 
These  factors  are  obtained  via  the  LINPACK  subroutine  SGECO.  As  a  byproduct  of  this 
subroutine  call,  an  estimate  of  the  conditions  numbers  of  Yj.  Ki  ^  /cCY*).  is  returned. 

Step  6  :  Compute  i  such  that  k-  =  min  {kJ. 

I  €  N 

Remarks  :  Since  Y;(k)  has  no  special  structure  and  an  estimate  of  the  condition  number  is  needed, 
the  LINPACK  subroutine  SGECO  is  the  obvious  choice  in  step  5.  Algorithm  EYEHAT  is 
computationally  intensive!  This  is  the  price  paid  for  basing  the  choice  of  which  matrix  to  invert  on 
the  lowest  condition  number  criterion.  Another  guideline  might  be  to  invert  the  matrix  with 
smallest  dimension.  This  would  lessen  the  computational  burden  born  by  the  algorithm,  but  might 
lead  to  inaccurate  results  in  unusually  conditioned  circumstances. 

It  is  clear  that  the  first  feedback  must  be  found  by  finding  an  inverse  of  an  arbitrary  matrix. 
Subsequently,  the  other  feedback  can  be  calculated  via  (3.2.10)  using  algorithm  MLTPLY.  The 
next  algorithm  computes  the  feedbacks  for  a  2-player  game  where  the  value  of  i  is  already 
available  from  the  criterion  defined  in  algorithm  EYEHAT. 

Algorithm  FEEDBK 

Step  1  :  Compute  T\—  |y;  W,  using  the  LINPACK  subroutine  SGESL  where  the  L-L 

decomposition  is  obtained  from  algorithm  EYEHAT. 

Step  2  :  Compute  U t  =  I  —  Yj(k )  K  j(k  +  l ).  j  ^  i  in  one  matrix  multiply  Do-Loop. 

Step  3  :  Compute  V  (  =  TlUl  using  algorithm  MLTPLY. 


.  S.  -  i  •  .  •  . 


Step  4:  Compute /^(k)  =  V-  A  using  algorithm  MLTPLY. 

Step  5  :  Form  A  j  =  A  —  B  •  F  ■  where  j  ^  i  in  one  matrix  multiply  DO-Loop. 

Step  6  :  Form  the  other  feedback  F  j  =  ^  A  j  where  j  ^  i. 

Remarks  :  Algorithm  FEEDBK  is  sufficient  for  computing  the  feedbacks  only.  However,  the 
ultimate  goal  is  to  perform  one  iteration  of  each  Riccati  equation  described  by  (3.2.8).  For  that 
reason,  it  is  desirable  to  modify  step  5  above  so  that  other  quantities  are  available.  Specifically, 
change  it  to  Form  A  x  =  A  —  B 2  F2  and  A2  =  A  B i  F i  in  one  matrix  multiply  DO-Loop. 


3.2.4  Computation  of  Coupled  Riccati  Iterations 

From  (3.2.8f),  notice  that  if  F,  —  0,  then  just  compute  the  Lyapunov  equation 

ET  K,(k)E  =  A[  (k)  A'jCk+l)  A  j(k)  +  Q{ 
using  algorithm  QDFORM.  Otherwise,  if  F i  =  0.  j  »«*  i  then  compute 


(3.2.15) 


ET  K,(k)  E  =  At  A-i(k+l)  -  %T  ^  %  A  +  (2j  (3.2.16) 

using  algorithm  DISRIC.  Regardless  of  the  condition  encountered,  the  following  algorithm 
performs  one  iteration  of  coupled  discrete-time  Riccati  equations. 


Algorithm  RJCCAT 


Step  1 :  Compute  the  feedbacks  using  algorithm  FEEDBK. 
Step  2  :  Iterate  each  Riccati  equation  using  algorithm  DISRIC. 


Remarks  :  Computation  of  Multirate  Coupled  Riccati  Iterations 

If  L  €  Z+  the  Riccati  equations  to  be  iterated  are  the  same.  However,  whenever  L  i  Z+.  several 


simplifications  are  possible.  First. 


F,{ k)  =  ¥j(k)  A 
F2( k)  =  0  . 

Hence,  the  Riccati  equations  that  have  to  be  iterated  reduce  to 


(3.2.17a) 


(3.2.17b) 


.»■  A  ■«*  ■h  A.Wl  A.’t  .U  .'i.i'*..  »■«  ‘  t'U'ti 


ET  K  t(k)  E  =  Ar  l^r^k+D-^/T!  ¥,  ]a  +  <2i  (3.2.18a) 

Et  K2(X)  E  =  A  2  (k)  A^2(k+1)  A  2(k)  +  Q2  .  (3.2.18b) 


Concluding  Remarks  : 

Section  3.2  describes  the  numerical  solution  of  a  coupled  discrete-time  Riccati  equation 
problem  motivated  in  Chapters  1  and  2.  During  the  development  of  the  solution  method,  several 
algorithms  are  defined.  Collectively,  they  serve  to  provide  a  rich  numerical  foundation  upon  which 
more  sophisticated  algorithms  can  be  built.  Thus,  the  Riccati  algorithms  presented  here  are  an 
example  of  using  software  to  efficiently  solve  a  frequently-formulated  game  problem.  This  is  the 
overall  spirit  behind  CACSD  endeavors.  The  results  discussed  in  this  subsection  apply  to  the 
standard  LQ  regulator  problem  and  the  Nash  equilibrium  solution  of  an  LQ  state-feedback  problem 
(the  coupled  Riccati  case).  The  issues  associated  with  solving  single  and  coupled  discrete-time 
Riccati  equations  are  delineated.  As  a  consequence  of  this  analysis,  it  is  shown  that  the  numerical 
intensiveness  is  directly  related  to  the  criterion  used  for  determining  i  via  algorithm  EYEHAT. 
Beyond  that  fact,  the  Riccati  iteration  procedure  basically  boils  down  to  computing  one  or  more 
quadratic  forms,  most  conveniently  calculated  by  algorithm  QD FORM. 


3.3  Finite  Horizon  Problems  —  Existence  Issues 


In  this  section  we  investigate  conditions  for  which  the  existence  of  solutions  to  finite-horizon 
problems  is  guaranteed.  To  begin  with.  Proposition  3.1  in  Section  3.1  is  based  upon  a  singular  value 
argument  applied  to  H,(k).  However  utilizing  the  definition  (3.1.7).  a  stronger  result  can  be  stated. 


Define  : 

n12(k)  £  ¥,(k)fl2 
n21(k)  t  'Mk)#!  . 

Then.  E^k)  =  I  —  fl12(k)  fi21(k)  and  E2(k)  =  I  —  fl21(k)  fll2(k) .  Hence. 


(3.3.1a) 


(3.3.1b) 


Proposition  3.2  :  Existence  of  Feedbacks  for  Finite-Time  Problems 

K  |(k)  exists  if  and  only  if  \(  fl12(k)  ft2i(k))  ^  1  for  all  eigenvalues.  X2(k)  exists  if  and  only  if 
A(  fl21(k)  fl12(k))  ^  1  for  all  eigenvalues. 


Proof  :  Apply  the  similarity  transform  that  puts  flj2(k)  021(k)  into  Jordan  form  to  sj(k)  and 
the  first  statement  follows  immediately.  Apply  the  similarity  transform  that  puts  fl2i(k)  ft12(k) 
into  Jordan  form  to  Hj(k)  and  the  second  statement  follows  immediately. 

□ 


3.4  Infinite  Horizon  Problems  —  Convergence  Issues 

Before  concluding  this  chapter,  an  investigation  of  the  existence  of  solutions  to  the  infinite¬ 
time  LQ  Nash  game  is  conducted.  To  the  author's  knowledge,  virtually  no  work  has  produced 
either  necessary  or  sufficient  conditions  regarding  existence  of  solutions  even  though  the  theory 
governing  the  one-player  optimal  infinite-time  LQ  regulator  problem  is  well-established  [5.48]. 
Solutions  to  single-rate  infinite-time  LQ  Nash  games  will  be  studied.  We  determine  conditions 
under  which  the  existence  of  an  infinite-time  solution  is  guaranteed. 


3.4.1  Preliminaries 

Let  X  denote  the  space  of  all  nxn  symmetric  matrices.  Let  Y  C  X  denote  the  set  of  all 
posit ive-semidefinite  symmetric  matrices  in  X  .  Clearly  X  is  a  linear  vector  space  and  Y  is  a  closed 
subset  of  X  .  Let  X  .  Y  €  X  be  any  two  arbitrary  elements  of  the  space  X  .  Then 


(X.Y)  t 


X  0 
O  >’ 


denotes  an  arbitrary  element  of  the  product  space  X  X  X  .  Now.  let  !  I2  denote  the  standard 
induced  matrix  2-norm.  Then  for  any  (X  .  )'  )  €  X  X  X  .  define 


II  (X  .  Y  )#  J  n  X  +  H  Y  I,  (3.4.1) 

to  be  the  norm  on  the  product  space.  Obviously,  the  space  X  X  X  with  the  norm  defined  by  (3.4.1  ) 


is  a  complete  normed  linear  vector  space.  Hence,  it  is  a  Banach  space. 


Assumption  3.1  :  Throughout  the  remainder  of  this  section,  the  matrix  E  is  always  assumed  to 
equal  the  identity. 

We  have  two  maps  : 

Jf,  :  XxX  -  X  and  R2  :  XxX  -  X 
Specifically,  from  (2.1.11c).  consider 

Rt  ( x.y )  *  A[(x.y)  [x^  +  fl,*,-1#?  |~1>i1(x.y)  +  Qi 


r2(x.y)  J  A£(x.y)  y-1 +52r2-1s£  ]  1  A2(x.y)  +  e2 

whereA,(X.y)  *  >1  -  fl2F2(X  .  y  )  and  A  2(X  .  Y  )  *  A  -  BxFx{X  .Y)  .  Define 


4>i(x)  *  x_1 +  5i^r1sf  1 . 


Then. 


If,  (X.y)  ^  A[(X.y)*,(X)A,(X.y)  +  <2, 
r2(x.y)  *  A£(x.y)<i>2(y)A2(x.y)  +  g2 


(3.4.2) 


Notice  that  A  x(X  .  Y  )  and  A2(X  .  y )  serve  to  couple  the  two  equations  of  (3.4.2).  Consider 
stacking  the  previous  two  equations.  Then  define  : 


R(x.y)  ± 


X, (x.y)  o 
o  x2 (x.y) 


Given  this  setup  we  observe  that 


R  :  XxX  -  XxX 


The  following  result  is  the  most  important  one  of  this  chapter.  It  provides  sufficient 
conditions  for  existence  of  and  convergence  to  the  single-rate  infinite-time  LQ  Nash  equilibrium 


solution. 


Theorem  3.1  ;  Existence  and  Convergence  of  Solutions  to  the  Infinite-Time  LQ  Nash  Game 
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Given  the  system  (2.1.1)— (2.1.2)  with  performance  indices  (2.1.3).  then  there  exist  constants. 
0  <  €7  <  1  and  0  <  €£  <  1  which  define  the  matrices  R  j  ^  I  and  R  2  ^  ~r  I  such  that  for  all 

R 1  and  R2  where  a(R  t)  >  -i-  and  ff(^2)  >  —  and  all  system  matrices  A  with 

€,  r2 

I A  I  <  £(  €7.  6^)  <  1.  the  set  of  coupled  difference  equations  : 

,  ,  -1 

K i(k)  =  ,4[(k)  (tfjOc+l))-1  +51(^,(k))*15{  A  ^k)  +  <2i 
A-2(k)  =  ^4 ^ (k)  U/k+l))-1  +  B2(R2{)/l)Y1Bt2  *  A2(k)  +  Q2 

constitutes  a  contraction  mapping  and  hence  converges  to  a  unique  fixed  point  denoted  ( K  \  .  K 2  )  . 


(3.4.3) 


Proof  :  To  begin  with,  let  €7  — *  0  and  — *  0  independently.  From  Theorem  2.1  we  know  that  in 
the  limit.  I Fi  I  =  I F2 1  =  0  which  implies  that  A  1  -*  A  and  A2~*  A  in  the  limit.  Furthermore. 
I(^’j)-1l  =  «7  -*  0  and  l(/T2)-**  =  -♦  0  by  construction.  Hence,  the  two  coupled  discrete-time 

Riccati  equations  tend  toward  two  decoupled  discrete-time  Lyapunov  equations  given  by 

*,(k)  =  At  *t(k+lM  +  <2i 

(3.4.4) 

K2{V)  =  AT  K 2(k+l )  A  +  Q2 
in  the  limit  as  €7  -*  0  and  -*  0  . 

It  is  well  known  that  this  set  of  equations  has  a  unique  positive-definite  symmetric  fixed  point 
(K  l  .  K 2  )  if  and  only  if  I  A  fl  <  1.  Moreover,  (3.4.4)  is  a  contraction  mapping,  so  that  for  any 
positive-semidefinite  symmetric  initial  guess  (K  j(0),  AN(O))  .  convergence  to  the  fixed  point  is 
assured  as  k  — *  —00. 


i 


1 


\ow  we  must  argue  that  for  a  given  problem,  there  exists  an  open  neighborhood  of 
R  =  R  2~^  =  0  characterized  by  the  constants  e7  and  €7  such  that  tor  all  A  with 
II .4  II  <  £(  €p  <  1.  equation  (3.4.3)  constitutes  a  contraction  mapping.  The  next  subsection 
will  more  than  fulfill  this  requirement. 

□ 


T' t;  V '  V 


H 


s 


f>: 


«.v 
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3.4.2  Contraction  Mapping  Argument 

If  the  mapping  R  is  to  be  a  contraction,  then  the  following  condition  must  be  satisfied  : 

I  R(X .  Y )  -  R(Z  .  W  )  I  <  a  I  (X .  Y )  -  (Z  .  W )  I  (3.4.5) 

for  all  X .  Y .  Z  .  W  €  Y  and  0  ^  a  <  1  .  Straightforwardly,  compute  : 

IR(X.  Y)  -  R(Z.  W)l  =  IX,  (X.K)  -  X,  (Z.W)l2  +  I  X2  (X.  Y )  -  X2  (Z  .  W)  ^ 
Similarly, 

l(X.K)  -  (Z.W)I  =  IX-ZIz  +  IY-W^. 

Therefore,  the  condition  (3.4.5)  is  equivalent  to 


ix,  (x.y)  -  x,  (z.w)i2  +  ix2  (x.n  -  X2  (Z.W)l2 
<  a|lX-ZI,  +  IF-WI, 


(3.4.6) 

The  purpose  of  the  ensuing  discussion  is  to  completely  characterize  those  cases  for  which  (3.4.6) 
holds.  The  contraction  mapping  that  we  seek  occurs  on  a  closed  subset  of  Y  X  Y  .  Hence  it  must  be 
shown  that  there  exists  a  region  A  C  Y  X  Y  where  (3.4.2)  maps  A  into  itself.  Toward  this  end.  the 
following  facts  are  established. 


Fact  3.1  :  Given  any  X  €  Y  .  then  I  4>,(X  )  l2  I  X  ^  for  all  i  €  {  1.2  }. 


Proof  :  Since  X  ^  0  and  XjXf^Xf  ^  0  . 


I  4>,(X  )  l2  = 


=  JX 


^(X-1  +  B,R-*B{  )  o;(X-1) 

where  cr(  )  denotes  the  smallest  singular  value  of  (  )  . 


To  simplify  the  derivation,  introduce  the  following  definitions 


p,(X.  }•  )  J  M,(X  .Y)U  ,  i=1.2 


8,(P,)  t 


»2 


i  -( p,(x.  n  )2 


□ 


(3.4.7) 


•■r 

*/ 


.  i=  1 .2 


( 3.4.S) 


54 


Bh  *  {X  eX  \  IX  l2^&)  .  (3.4.9) 

Equation  (3.4.9)  describes  a  closed  ball  of  radius  8  . 

Then  we  have  : 

Claim  3.1  ;  Existence  of  a  Region  where  the  Riccati  maps  are  Into 

Suppose  p;(X .  Y  )  ^  J>x  <  1  and  f  - — '  —  —  for  i=l  .2  .  Given  any  (X .  Y )  €  Y  X  Y  such 

1  -  (  p  i  )2 

that  I  X  k  ^  and  I  Y  §2  ^  T2  .  then  for  all  such  X .  Y  it  follows  that  (X ,  Y  )  ^  R  (X .  Y  )  has  the 
property  that  I  X  l2  ^  and  I  Y  ^  ^  ^2  •  Hence.  R  (X .  Y  )  maps  A  ^  B-gt  X  Bj2  C  Y  X  Y  into 
itself. 


Proof  :  Notice  that  from  (3.4.2)  : 

ix  i2  =  iAr1(x.Y)*1(x)Al(x.Y)  +  q,i2 

<  lAxiX.Y)#  I4>1(x)l2  +  lfi,l2 

^  iAjX.Y)!?  IX  Ij  +  I(2i«2 
where  Fact  3.1  has  been  utilized  in  the  last  step. 

Thus. 


ix  «2  «;  (P,(x.y))2  ix  i2  +  i<2, i2 


^  (p,)2  Si  +  lQll2 


Similarly. 


(Pl)2»gl»2 

1  ~  (  Pi  )2 


+  «Gl«2 


=  5,. 


I  Y  ll2  =  WA^(X.Y)  4 >2(y  )A2{X,Y)  +  Q2  ll2 

^  ».4,(x.  n»22  it  i2  +  »<22»2 


^  (p,)2  8,  +  #Q2« 2 
(p  2)2  II  (2  2  «2 


I  -  ( p ,  )2 


+  n  <2  2  h2  =  "F,  . 


Therefore,  I  X  »3  ^  8,  .  II  Y  h  <  53  .  and  (X .  Y  )  €  A  f  Bjt  X  Bz2 


I 


□ 


i 
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I 

8 


Now.  we  proceed  to  derive  conditions  such  that  (3.4.6)  holds  in  the  region  A  .  In  particular, 
focus  on  the  map  Rt  (X .  Y) .  A  parallel  argument  can  be  developed  for  the  map  R2  (X.Y). 
Given  any  (X.Y)  and  (Z  .  W )  €  X  X  X  , 


Ri  ( X.Y )  -  Rt  ( Z.W ) 

=  Arl(X.Y)<t>1(X)Al(X.Y )  -  Af(Z.  W  )  4>,(Z  )  A  ,(Z  .  W) 
Application  of  the  matrix  identity 


(3.4.10) 


£ 

"i. 


S. 

'w 


9! 


v. 

v. 

V 


W  X  -  Y  Z  =  -L 
2 

to  (3.4.10)  yields : 


(W  -  r)(X  +Z)  +  (W  +  Y)(X  -Z) 


(3.4.11) 


Rx  (X.Y)  -  Rt  (Z.W) 


2 


|A[(x\n-A[(z.vn|-  [<t,(A:)A1(x.r)  +  <i>1(z)A1(z. w) 


+  \a[(X.Y)  +  At1(Z.W)  ]  •  |  <t>x(X )  A  i(X .  Y  )  -  4>X(Z  )  A  ,(Z  .  W )  | 


\b2(F2(Z.'W)-F2(X,Y)  )  jr  •  [<t>i(X)Al(X.Y)  +  4>l(Z)A1(Z.W)  j 


2 


+  |  A\(X.Y)  +  A\(Z.W)  |  |<J>1(X)A,(X.r)-<J>1(Z)A1(Z.  VV) 

But.  from  an  additional  use  of  (3.4.11) 

*l(X)Al(X.Y)-*l(Z)Al(Z,W) 


(  4>,(X  )  -  <D,(Z  )  j  (a1(X.K)  +  A,(Z.W)  | 


+  *,(X)  +  fl>1(Z)  |-  \AtiX  .Y)-Ai(Z.W) 

Therefore,  we  conclude  that  : 


/. 

A 


i 


(X.Y)  -  Rt  (Z.W) 

[  B2(  F  2\Z  .  W)-  F2(X  .Y)  ) 


<t>,(X  )  A  ,(X  .  Y)  +  <J),(Z  )  A  ,(Z  .  W  ) 


S 
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+  l^(x.n^Kz.'v))  |ti(x)_,i(z)||y(iCx.n  +  ,4i(z.w)) 

+  [‘4‘U'r)^'<Z>t')  1  ■  [*,(X)  +  »,(Z>  |  |  it2  (f2(Z,  HO -fi(X.n)  | 

Applying  the  matrix  2-norm  and  then  the  triangle  inequality  on  the  last  equation  gives 


I  j?i  (X.  K)  -  R  i  (Z.  W)  I2 

<  [ } A l(x ; Y )  *  'JOLUl  •  1  b2 ( f2(z .  w)- /-2(x . y ) ) *2  + 
"  2 

MrfX.n*a,(Z.lir)|?.,fcrif,_^,L  + 

I A  t(X  ,Y)  +  A  j(Z  .  IV )  12  ( j  +  ^^2  )  ^  .  |  Bl  (  ^2(Z  ,  W )  -  /^(X .  T )  )  «2 


(3.4.12a) 


The  corresponding  inequality  for  the  map  J?2  (X  .  Y  )  is 
1X2  (X.X)  -  X2  (Z.  W)J2 

I  4>2(K  )  A  2(X  .  y  )  +  <1>2(  w  )  A  2(Z  .  W  )h_  ■!  B  (  p  (7  vp)-F,(X.r)  )|7  + 

^  2 

M,(x.n  +  xa(z.w)i;2 , aJwr>u+  (3412b) 

4 

m2(x.  y )  +  a2(z  ,  w  )n2  fl }  +  ^2(vv } ,2  1 5, ( f,(z . w ) - f,(x . y ) ) i2 

4 

In  order  to  produce  the  conditions  for  which  (3.4.6)  holds  in  the  region  A  _  5?,  X  Z*52  . 


R,  J  JL  I  .  €,  >  0  .  i  €  {1.2} 


(3.4.13) 


and  show  that  for  any  and  all  positive-semidefinite  symmetric  X,  with  giR ,)  ^  ^ y  >  0 


II  R,(X.  y  )  -  /?,(Z  .  W)  ll2  <  a(  8.)  •  I »  X  -  Z  l2  +  II  y  -  W  «2  (3.4.14) 

where  (X  .  Y  ).  (Z  .  W )  €  A  .  0  <  a  (  S',)  <  1  .  and  T  (  5,)  is  yet  to  be  determined.  Then,  the 

contraction  mapping  argument  will  become  obvious. 

Therefore,  some  useful  facts  are  staled.  The  proofs  may  be  found  in  Appendix  B.  The 
following  result  introduces  the  important  variables  vx  and  v2  ■ 


v  .  v v  s.*;  v' 


Fact  3 J2  :  Given  any  X .  Y  6  and  R{  £  —  I  .  €i>0.  then  whenever  €t  <  - - -  holds  for 

'  €i  ^ ,  »22 

all  i  6  {1.2},  I  4\(X )  -  ^(r )  ^  ^  - 1 - I  X  -  Y  ^  where  v,  *  e,  8*  I  B,  l?  <  1  . 

(1-Vi)2 

Proof  :  See  Appendix  B 


Fact  3.3  :  Given  any  (X.  Y)  €  A  ^  Bix  X  Bjt  and  Ry  ^  —  I .  €*>0,  then  whenever 
0  <  €j  <  ej  (^)  ^ T  holds  for  all  i  €  {1.2}.  it  follows  that 

p,(x . r )  =  i a j(x . y ) h  ^  t  ■  *i  >  1 

p2ix.r )  =  ia2(x,y)i2  ^  £  1  .  >ci  ^  i 

1  —  Vi  v2 

where  ^  *  e^,  S  Bt  \}  and  £  *  II  A  l2  . 


Proof  :  See  Appendix  B 


Consider  (3.4.12a)  and  (3.4.12b)  where  (X  .  Y  ).  (2  ,  W  )  €  A  .  Consequently.  (3.4.12a)  coupled 
with  Facts  3.1.  3.2.  and  3.3  yields  : 

i*i  (x.r)  -  *,  (z.w) n2 

<  1  *  'YL ■  'l  1± ‘(z : > ' !i  ■  ■  b,  (  f  ,(z .  w  >  -  F  ;(x .  n  > + 

-— 2 - *'  IX-ZI,  + 

4  (  1  —  )2 


2£  Ky 


X  ll2  +  II  Z  »2  I  1  B2  (  F2{Z  .  W  )  -  F2(X  .  Y  )  )  ll2 

(  £  »T  I2 

n  n  /¥-•/  '7  iir  \  rT/vx,\'\n  i  I  *  1  ) 


<  £  irx  ii  b2  (  f2(z  .w)-  f2(x  .  r ) )  i2  +  —  L  -ix  -z  u  + 

1  Uy 

£8y  Ky  I  B2(  F2U  ,W)~  F2(X  .  Y  )  )  H2 


E 


S 
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^  2^*1  I52(F2(Z,  W)-F2(X,y))i2  + 


£  «i 

1  -  V~l 


2 

I  x  -  Z  l2 


(3.4.15) 


and  (3.4.12b)  becomes 


I  if  2  ( X.Y )  -  R2  (Z.W) 

2fJ2ici  I Bl{Fx(Z.W)-Fl{X.Y'})%2  + 


£  *2 
1  -  v\ 


2 

I  Y  -  W  l2 


(3.4.16) 


Equations  (3.4.15)-(3.4.16)  are  nearing  the  form  of  (3.4.14).  The  next  result  allows  the 
contraction  mapping  argument  to  be  completed.  The  details  of  the  proof  may  be  found  in 
Appendix  B. 

Theorem  3.2  :  Lipschitz  Constants  for  the  Feedbacks 
Given  (X.Y).  (Z  .  W)  €  A  ,  then 


F1{Z.W)-F}(X.Y))  Ij  <  at1}  1  X  —  Z  l2  +  <*,2  I  T  ~  W  l2 

ii  s2  ( f2(z .  w )  -  f2(x  .  y ) )  i2  <  a2i  i x  -  z  n2  +  <*22 1  r  -  w  n2 

£v\k\k2  £  *1^2*1  *~2  £v\V2K\*2  .  £  *2  *~1  *2 

where  au  =  - - -  .  a12  =  - - -  .  ar21  =  - - -  .  and  a22  =  - - - 


Proof  :  See  Appendix  B 


□ 


Making  use  of  Theorem  3.2.  (3.4.15)  may  be  rewritten  as  : 
lie,  (X.Y)  -  R}  (Z.  W)ll, 


2^8,^ 

a,,  II  X 

~  7  ll2  + 

<*\ 

2  |  5,  k] 

a,,  + 

£  ><\ 

2 

i  -  irx 

£  x\ 


X  -  7.  IU 


1  -  V  i 

X  -  Z  «2  +  2  £1$,  k]  al2  II  >'  -  W 


1 

) 


I 
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=  2  i2  vx  (  K\)2  k2  +  **  (  *l] *  ■  I  X  —  Z  I2  +  2  i2  Ji  r^OTi)2^  ir  -W  Iz 

1  -FT2  ^2 


=  I2  (  Kp2  liTiiCj  +  - 1 - _  I  A  -  Z  I2 

(l-vl)2 


T 

C2  ,-7-  ,T  (  =-\2  s- 


+  2  £2  ^(Kp^  iy-wij 

S2 


Likewise.  (3.4.16)  reduces  to 


1*2  (x.n  -  *2  (Z.  W)l2 


<  2  £2  — i  Vy^KyU^)2  XX  ~Z  ^ 
51 


+  £2  (  *T)2  2  ^  /q  +  - ! - _  -|  Y  -  W  #2 

(  1  -  ^  )2 


Finally,  adding  (3.4.17a)  and  (3.4.17b)  yields  the  desired  result. 


*,  (X.Y)  -  *,  (Z.W)U  +  1*2  (X.l’)  -  *2  (Z.  W)ll2 


(  — ^2  "T" 

^  f2  2  FT  (  kJ)2  kT  +  - 1 - -  +  2  -3-  F]  F]i  *T  (  K~2y  ■  0  x  —  Z  <2  + 

(  1  -  FT  )2  T, 


I2  2  Fi  iq  (  icT)2  +  — — .  +2  Ji  FT  FT  (  Kp2  *;  II  y  -  W  I2 

(  1  -  f;  )2  ^2 


I!  ?  ax  ■  X  X  -  Z  X2  +  £2  c*2  •  II  Y  -  W  1, 


^  f-  ■  max  {a,.  a2)  II  X  -  Z  I,  +  II  Y  -  W 


Clearly.  (3.4.18)  and  (3.4.6)  are  of  the  same  form  with 


a  =  £2  max  {<*1.  a2| 


(3.4.17a) 


(3.4.17b) 


(3.4.18) 


(3.4.19) 


Moreover.  £  1  A  ll2  can  always  be  chosen  sufficiently  small  enough  so  that  0  ^  a  <  1  .  Hence,  the 


existence  of  a  non-triviai  region  where  coupled  discrete-time  Riccati  equations  behave  as  a 


contraction  mapping  is  established. 


It  is  apparent  that  there  are  3  parameters  which  govern  the  region  where  a  contraction  mapping 


»  -  •  *  4  *  .«  *4  ~  * 

‘V.V '>> '  •  V- 
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occurs.  These  are  ,  €2  •  and  £  Before  concluding  this  section,  a  procedure  is  developed  for 
determining  those  ranges  of  values  for  which  a  contraction  mapping  is  guaranteed.  Essentially, 
these  quantities  must  be  chosen  small  enough  so  that  all  the  assumptions  of  the  derivation  remain 
valid.  To  begin  with,  £  is  required  to  be  smaller  than  some  £  as  yet  to  be  determined.  From 
(3.4.7)  and  Claim  3.1  : 


Pi(X.T)  *  lAjLX.Y)^  *  W;  <  1  • 

But  in  view  of  Fact  3.3  and  without  loss  of  generality  : 


(3.4.20) 


Pi  t  £*>  <  1 


(3.4.21) 


provided  £  <  min  -L,  -ir-  and  0  <  €.  <  e,  ( 1>.)  ^  - ? — —  .  i  =  1,  2.  However,  there  is  a  more 

*i  *2  “  I,lfi,»22 

restrictive  condition  imposed  on  £  by  (3.4.19).  That  is 


£  <  mini 


_J _ 1_ 

V<* 1  V« 2 


(3.4.22) 


1  1 


In  fact,  it  is  easilv  shown  that  min  - .  -  <  min  — .  —  with  equalitv  if  and  only  if 

Val  lK‘  K2 


=  €■>  =  0  .  Therefore,  define  £  ^  min  — - — .  — 1 — 


and  pick  £  such  that  £  <  £  . 


Now.  a  moment's  reflection  reveals  that,  in  general,  the  3  parameters  €j  .  e2  .  and  £  cannot  be 
explicitly  solved  for.  To  see  this,  consider  £  which  must  be  chosen  smaller  than  £  .  Well.  £  is  a 
function  of  a!  and  a2  which,  in  turn,  are  functions  of  Sj  and  b2  through  (3.4.18)  .  But.  "5";  is  a 
function  of  p,  by  (3.4.8)  and  (3.4.21)  indicates  that  p",  is  a  function  of  £  .  Hence  £  is  a  nonlinear 
function  of  itself.  Similar  arguments  can  be  stated  for  €j  and  e2  . 


Therefore,  the  task  at  hand  is  to  enumerate  the  cases  where  6)  ,  e2  ,  and  £  yield  a  contraction 
mapping.  The  following  algorithm  accomplishes  this  task. 


Contraction  Mapping  Surface  Generator  Algorithm 
Step  1  :  f  or  v\  .  €  (0,  I  )  but  fixed  do  the  following  : 
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i 


! 


:■£ 


& 


$ 


Step  2  :  Compute  Kj  and  k2  as  defined  in  Fact  3.3. 
Step  3  :  Choose  £  =  min 


_1_  J_ 
*1  ’  *2 


Step  4 :  For  £  6  (0,  £),  compute^  (£)  and^2  (£)  • 
Step  5  :  Compute  and  a2  as  defined  in  (3.4.18)  . 


Step  6  :  If  £  <  min 


_1 _ 1_ 

Vo i  -v/a 2 


then  pick  €2  and  e2  to  give  and  v2  as  fixed  in  Step  1. 


Save  the  triple  (€j.  e2,  £)  as  a  valid  combination. 
Step  7  %  If  not  done  then  go  to  Step  1  else  stop. 


A 


S" 

.v 

V 


a 


pi 


u 


If  this  algorithm  is  implemented  on  a  computer,  then  a  list  of  possible  maximal  values  can  be 
compiled.  Actually,  the  roles  of  v]  and  v~2  are  completely  interchangeable  insofar  as  the  norms  of 
Q 1  and  Q 2  are  equal.  Table  3.1  summarizes  various  maximal  ranges  of  the  3  parameters. 

For  example,  let  I  B  ,  l2  =  U?2*2  =  *Qih  ~  =  1  Given  any  and  R2  with 

a(tf ,)  .  a(R2)  >  q  4^qq  =  2.0833  and  any  A  with  £  =  W(A  )  <  0.1  ,  then  the  corresponding 

coupled  discrete-time  Riccati  equation  iterations  will  be  a  contraction  mapping  for  all  initial 
conditions  (  X„.  F0)  where  cr(X0)  .  <f(K0)  ^  1.04  .  Note  that  if  dynamic  programming  were  used 
to  solve  these  equations,  then  the  initial  condition  would  be  (  X<>.  F0)  =  (  Q  i-  Qt)  and  hence  the 
requirement  that  the  initial  conditions  lie  in  a  ball  of  radius  1.04  times  the  magnitude  of  (2i  is 
automatically  satisfied.  Furthermore,  the  contraction  mapping  constant  defined  by  (3.4.19)  is 
a  <  (O.l)^  28.0  =  0.28  .  Obviously,  the  requirement  that  cr(A  )  <  0.1  is  a  conservative  one. 
This  is  due.  in  part,  to  the  large  discretization  increment  of  Table  3.1. 


(V 

i 


'  « 


3.4.3  A  Numerical  Example 

In  this  subsection,  we  construct  a  numerical  example  that  illustrates  the  contraction  mapping 
derived  in  the  last  subsection.  Consider  the  system  (2.1.1 )-( 2.1.2)  with  performance  indices 
(2.1.3).  Let  Ct  =  C 2  —  I-  =  -S  [  =  S2  =  I  Then  choose 


Table  3.1.  Normalized  Maximal  Ranges  of  Values  for  Which  the  Coupled  Riccatis  Are  Contractions 


*2 

£ 

81 

S2 

€1 

€2 

<*1 

<*2 

0.10 

0.10 

0.70 

2.53e+00 

2.53e+00 

0.0395 

0.0395 

1.83e+00 

1.83e+00 

0.10 

0.20 

2.17e+00 

1.83e+00 

0.0460 

0.1093 

2.24e+00 

2.67e+00 

0.10 

0.30 

0.50 

1.82e+00 

1.47e+00 

0.0551 

0.2035 

2.71e+00 

3.81e+00 

0.10 

0.40 

0.40 

1.52e+00 

1.27e+00 

0.0660 

0.3160 

3.24e+00 

5.4U+00 

0.10 

0.50 

0.30 

1.29e+00 

1.14e+00 

0.0776 

3.84e+00 

7.81e+00 

0.10 

0.60 

0.20 

1.13e+00 

1.06e+00 

0.0884 

4.52e+00 

1.18e+01 

0.10 

0.70 

0.20 

1.15e+00 

1.06e+00 

0.0866 

0.6608 

5.24e+00 

1.97e+01 

0.10 

0.80 

0.10 

1.04e+00 

1.01e+00 

0.0962 

0.7886 

6.08e+00 

4.10e+01 

0.20 

0.20 

0.50 

1 .64* +00 

1.64e+00 

0.1219 

0.1219 

3.38e+00 

3.38e+00 

0.30 

0.40 

1.44e+00 

1.35e+00 

0.1388 

0.2218 

4.22e+00 

4.99e+00 

■KwaSl 

0.40 

0.30 

1.26e+00 

1.18e+00 

0.1583 

0.3388 

5.21e+00 

7.31e+00 

0.20 

0.50 

0.30 

1.33e+00 

1.19e+00 

0.1500 

0.4200 

6.35e+00 

1.09e+01 

0.60 

0.20 

1.15e+00 

1.08e+00 

0.1736 

0.5554 

7.73e+00 

1.68e+01 

0.20 

0.70 

0.10 

1 .04e+00 

1 .02e+00 

0.1922 

0.6864 

9.34e+00 

2.86e+01 

0.20 

0.80 

0.10 

1 .05e+00 

1.02e+00 

0.1908 

0.7837 

1.12e+01 

6.02e+01 

0.30 

0.30 

0.30 

1.23e+00 

1 .23e+00 

0.2449 

0.2449 

6.44e+00 

6.44e+00 

0.30 

0.40 

0.30 

1.29e+00 

1.24e+00 

0.2317 

0.3214 

8.21e+00 

9.77e+00 

0.30 

0.50 

0.20 

1.14e+00 

1.10e+00 

0.2626 

0.4532 

1.04e+01 

1.50e+01 

0.30 

0.60 

0.20 

1.18e+00 

l.lle+00 

0.2543 

0.5397 

1.31e+01 

2.39e+01 

0.30 

0.70 

0.10 

1.05e+00 

1.03e+00 

0.2861 

0.6810 

1.64e-t-01 

4.15e+01 

0.30 

0.80 

0.10 

1.06e+00 

1  -03e+00 

0.2832 

0.7766 

2.04e+01 

8.90e+01 

0.40 

0.40 

0.20 

1.13e+00 

l.l3e+00 

0.3556 

0.3556 

1 .29e+01 

1.29e+01 

0.40 

0.50 

0.20 

1.16e+00 

1.14e+00 

0.3438 

0.4388 

1.69e+01 

2.05e+01 

0.40 

0.60 

0.10 

1.05e+00 

1.04e+00 

0.3823 

0.5796 

2.22e+01 

3.37e+01 

0.40 

0.70 

0.10 

1.06e+00 

1.04e+00 

0.3777 

0.6735 

2.91e+01 

6.07e+01 

ran 

0.50 

0.10 

1.04e+00 

1.04e+00 

0.4800 

0.4800 

2.80e+01 

2.80e+01 

0.60 

0.10 

1.06e+00 

1 .05e+00 

0.4739 

0.5724 

3.83e+01 

4.81e+01 

0.50 

0.70 

0.10 

1.07e+00 

1.06e+00 

0.4658 

0.6627 

5.27e+01 

8.99e+01 

0.60 

0.60 

0.10 

1.07e+00 

1 .07e+00 

0.5625 

0.5625 

6.91e+01 

6.91  e+01 

*  Assumes  1  2 1  l2  =  0(2  2^2  =  *  =  I  5 1 12  =  $  B  2  ll2 


0.4755  0.0459  -1.13e-4 

A  =  0.0459  0.345  -7.175e-5 

— 1.13e  —4  -7.175* -5  0.25 


9.87e  2  1.37 

=  1.23  and  =  -1.0e-3  . 

-1.01e-3  ‘  1.0e-5 

The  eigenvalues  of  A  are  0.49.  0.33,  and  0.25.  Next,  let  R  t  =  R  2  =  11.0.  From  Table  3.1  it  is 
observed  that  this  case  falls  within  the  region  of  contraction  mapping.  To  illustrate  the  contraction 
mapping  behavior,  this  problem  is  run  for  14  iterations  us  ng  the  L-A-S  operators  described  in 
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Chapter  4.  A  full  listing  of  this  run  may  be  found  in  Appendix  C.  The  important  results  are 
summarized  in  Table  3.2.  The  values  of  A ,  and  A  2  at  the  end  of  the  iterations  are 

1.00233  1.76897e  —2  -3.43357e-5 

Ki  =  1.76897* -2  1.13453  -5.14615* -5 

-3.43357*  -5  -5.14615*  -5  1.06667 

1.00233  1.76895e  —2  -3.43351 e-5 

*2  =  1.76895* -2  1.13453  -5.14618* -5  . 

-3.43351* -5 -5.14618* -5  1.06667 

Observe  that  the  iterations  converge  rapidly.  After  14  iterations.  Aj(k)  and  A2(k)  are 

changing  by  no  more  than  1.0e-13.  Even  more  remarkable  is  the  fact  that  a  is  almost  constant 

throughout  the  iterations.  One  final  thing  to  note  is  that  the  predicted  contraction  mapping 

2 

constant  is  (0.5)  •  3.38  =  0.845.  while  the  observed  constant  is  almost  an  order  of  magnitude  less 
than  that.  In  general,  this  large  difference  occurs  because  the  inequalities  used  in  proving  the 
contraction  mapping  are  not  all  tight  simultaneously. 

Since  the  contraction  mapping  argument  guarantees  convergence  to  a  unique  fixed  point  and 
the  values  of  A-!  and  A2  are  accurate  to  1.0e-13.  then  we  conclude  that  these  values  for  K i  and  K2 

Table  3.2.  Results  of  Contraction  Mapping  Iterations 


Iteration  # 

AA  i 

aa2 

a 

0 

1.00000e+00 

1 .00000e-00 

1.20666e-01 

1 

1.20666e-01 

1.20666e-01 

1.18329e-01 

2 

1.42782e-02 

1.42782e-02 

1.18053e-01 

3 

1.68558e-03 

1.68558e-03 

1.18020e-01 

4 

1.98932e-04 

1.98933e-04 

1.18016e-01 

5 

2.34772e-05 

2.34773e-05 

1.18016e-01 

6 

2.77068e-06 

2.77069e-06 

1.18016e-01 

7 

3.26983e-07 

3.26986e-07 

1.18016e-01 

8 

3.85891e-08 

3.85894e-08 

1.18016e-01 

9 

4.554 12e-09 

4.55416e-09 

1.18016e-01 

10 

5.37457e-10 

5.37463e-10 

1.18016e-01 

11 

6.34284e-l  1 

6.34291e-l  1 

1.18016e-01 

12 

7.48554e-12 

7.48565e-12 

1.18014e-01 

13 

8.83404e-13 

8.83404e-l  3 

1 . 18020e-01 

14 

1.04246e-13 

1.04273e-13 

1.18080e-01 

:wy>> . 
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satisfy  the  Algebraic  Riccati  Equations  (ARE)  to  13  significant  digits.  Thus,  for  some  cases,  the 
coupled  Riccati  iteration  algorithm  provides  a  method  for  obtaining  the  solution  to  two  coupled 
discrete-time  AREs.  This  solution  may  also  be  the  solution  to  the  infinite-horizon  LQ  Nash  game 
problem. 


CHAPTER  4 


L-A-S  OPERATORS  FOR  SINGLE  AND  COUPLED  DISCRETE-TIME  RICCATI  ITERATIONS 

This  chapter  describes  the  L-A-S  [32-36]  operators  created  for  the  task  of  iterating  single  and 
coupled  discrete-time  Riccati  equations.  There  are  a  total  of  six  new  operators  to  the  L-A-S 
package.  As  presented  in  this  chapter,  they  are  SYST,  LQ,  DRE.  GAME,  LQNG,  and  MLTR.  Single 
Riccati  iterations  are  addressed  first.  Then  the  coupled  (game)  case  is  described.  Operators  SYST 
and  DRE  fall  into  the  first  category,  while  operators  GAME,  LQNG.  and  MLTR  fall  into  the  second 
category.  L-A-S  operator  LQ  is  used  for  both  single  and  coupled  Riccati  iterations. 

First,  a  brief  description  of  the  operator  is  given.  Next,  the  corresponding  excerpt  from  the 
L-A-S  Help-File  is  presented.  Then  a  typical  example  of  the  usage  of  the  operator  is  demonstrated. 

4.1  L-A-S  Operator  SYST 

The  L-A-S  operator  named  SYST  is  used  to  define  a  linear  shift- invariant  descriptor  system. 
Often,  it  is  the  first  operator  issued  to  the  L-A-S  interpreter  when  a  linear  quadratic  regulator 
problem  is  being  studied. 

L  —A  —S  Help  —File  Description. 

SYST  -  Descriptor  SYSTem  description 
Syntax  :  A.  B.  C  [.  D  [.  E]  ]  (SYST)  = 

Input  Data  :  A  [N.N]  ,  B  [N.P]  .  C  [M.N]  , 

D  [M.P]  .  E  [N.N]  (these  last  two  arrays  optional) 

Options  :  E,  L.  T 

Description  :  Identifies  the  following  discrete-time  descriptor  system 

E  .t(k+l)  =  .4  jc  ( k. )  +  B  u(k) 

>(k)  —  C  x(k)  +  D  u(k) 


Note  :  If  E  is  omitted,  it  is  assumed  to  be  the  identity  matrix. 
If  D  is  omitted,  it  is  assumed  to  be  a  zero  matrix. 
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Example  of  Usage 

An  example  of  the  usage  of  operator  SYST  is  given  below.  The  following  segment  is  an  L-A-S 

program  that  identifies  a  discrete-time  descriptor  system. 

:  Descriptor  System  Identification 

(inp)=a,b.c.d.e 

a.b.c.d.e(out)» 

a.b.c.d.e(syst)** 

Notice  that  the  program  is  completely  generic  in  the  sense  that  the  dimensions  of  the  matrices  a.  b. 
c.  d.  and  e  are  not  specified.  This  is  an  important  feature  of  the  L-A-S  language.  The  same  program 
can  be  run  over  and  over  again  using  different  matrices  (of  different  order)  each  time.  If  this 
program  is  executed  for  a  simple  second-order  system,  the  following  output  would  result. 


> :  Descriptor  System  Identification 
>(inp)=a.b.c.d,e 
***  Matrix  a  *** 

Enter  the  dimensions  of  this  matrix.  >2.2 

Matrix  :  a  Enter  C.D.E.I.N.P.R.Z  or  H  for  Help.  >R 
ROW  1  >1.2 
ROW  2  >3.4 


***  Matrix  b  *** 

Enter  the  dimensions  of  this  matrix.  >2.1 

Matrix  :  b  Enter  C.D.E.I.N.P.R.Z  or  H  for  Help.  >C 
COL  1  >1.0 

***  Matrix  c  *** 

Enter  the  dimensions  of  this  matrix.  >1.2 

Matrix  :  c  Enter  C.D.E.I.N.P.R.Z  or  H  for  Help.  >R 
ROW  1  >0.1 

***  Matrix  d  *** 

Enter  the  dimensions  of  this  matrix.  >1.1 
Enter  the  scalar  :  d  >0 
***  Matrix  e  *** 

Enter  the  dimensions  of  this  matrix.  >2.2 


Matrix  :  e 


Enter  C.D.E.I.N.P.R.Z  or  H  for  Help.  >1 


K 
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>a,b,c,d,e(out)» 


t; 

I 


1.000  2.000 

3.000  4.000 

b 

1.000 


i 

0. 

c 

0. 

1.000 

c*. 

<4 

d 

0. 

e 

1.000  0. 

0. 

1.000 

I 

>a,b.c.d.e(syst)’ 

The  matrices  d  and  e  are  optional  as  indicated  in  the  Help-File  description.  However,  if  matrix  e  is 
specified  then  matrix  d  must  also  be  provided  to  serve  as  a  placeholder.  That  is.  a  statement  such 
as 


A 

ft 

•I 


•> 

*•% 


a,b,c,.e(syst)= 

is  not  permitted. 


4.2  L— A—S  Operator  LQ 

The  L-A-S  operator  named  LQ  is  used  to  define  the  weighting  matrices  in  a  discrete-time 


3 

.  J 


system  or  game  problem.  Often,  it  is  the  second  operator  issued  to  the  L-A-S  interpreter  when  a 
linear  quadratic  regulator  or  Nash  game  problem  is  being  studied. 


L  —A  —S  Help  —File  Description. 


v 


LQ  -  Linear  Quadratic  weighting  matrices  for  system  or  game  theoretic  problems 

( >  r-\  «  n  4  r  /A  n  A  1  /  *  \ 


Syntax 
Input  Data 


Options 

Description 


Ql.  R1  [.  Q2.  R2]  (LQ)  = 

Q1  [N.N]  .  Rl  [Pl.Pl]  .  Q2  [N..N]  ,  R2  [P2.P2] 

(These  last  two  arravs  only  required  for  game) 

E.  L.T 

Identifies  the  weighting  matrices  in  a  linear-quadratic  discrete-time 
system  or  game  problem 
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Example  of  Usage 

An  example  of  the  usage  of  operator  LQ  is  given  below.  Consider  the  discrete-time  system 
defined  in  Section  4. 1 .  The  following  segment  is  an  L-A-S  program  that  identifies  the  same  system 
and  then  defines  two  weighting  matrices  in  preparation  for  the  study  of  a  linear  quadratic  regulator 
problem.  Rather  than  inputting  the  matrices  from  the  keyboard,  the  RDF  (Read  Data  File) 
operator  is  used. 


;  Linear  Quadratic  Regulator  Problem 
:  System  and  Weighting  Matrix  Identification 

:  System  Definition 
(rdf)-a.b.c 
a.b.c(out)« 
a.b.c(syst)- 

:  Weighting  Matrices  Definition 

(rdf)=q.r 

q.r(out)” 

q.r(lq)- 

If  this  program  is  run  using  the  second-order  system  of  Section  4.1.  the  following  output  results. 


> :  Linear  Quadratic  Regulator  Problem 
> :  System  and  Weighting  Matrix  Identification 
> ; 

> ;  System  Definition 
>(rdf)=«a.b.c 

Enter  name  of  the  Data  File  (DF)  for  matrix  a  >syst 

Opening  file  named  :  syst.DF 

Reading  array  named  :  a 

Reading  array  named  :  b 

Reading  array  named  :  c 

>a.b.c(out)= 

a 

l.(KK)  2.000 
3.000  4.000 

b 

1.000 

0. 


CJ. 


•y 


«  i 


y, 

■r. 


:  *, 


v. 


K 

Sc 


r 
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0.  1.000 


>a.b.c(syst)= 


> :  Weighting  Matrices  Definition 
>(rdf)«q,r 

Enter  name  of  the  Data  File  (DF)  for  matrix  q  >  lq 
Opening  file  named  :  lq.DF 
Reading  array  named  :  q 
Reading  array  named  :  r 

>q.r(out)- 

q 

1.000  0. 

0.  2.000 


1.000 


>q.r(lq)= 


4.3  L— A— S  Operator  DRE 

The  L-A-S  operator  named  DRE  is  used  to  iterate  a  single  discrete-time  Riccati  equation. 


L  —A  —S  Help  —File  Description 
DRE  -  Discrete-time  Riccati  Equation  iteration 


Syntax 
Input  Data 


K  (DRE)  =  KNEW 
K  [N.N] 


Output  Data  :  KNEW  [N.N] 

Options  :  E.  L.  T 

Description  :  Assuming  that  operators  SYST  and  LQ  have  been  issued. 

this  operator  performs  one  iteration  of  the  discrete-time 
Riccati  equation  defined  by 

Er  KNEW  E  =  Ar(K-Vrr'lOA  +  Q 

where  T  =  R  +  Br  K  B  .  and  ¥  =  (D"1  Br  K . 

A  ,  B  .  and  E  are  identified  by  the  SYST  operator.  Q  and  R  are 
identified  by  the  LQ  operator. 


-vk'lvf 


WHS 


sw 


WTO 


«:v 
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Example  of  Usage 

An  example  of  the  usage  of  operator  DRE  is  given  below.  Consider  the  discrete-time  system 
defined  in  Section  4.1.  The  following  segment  is  an  L-A-S  program  that  identifies  the  same  system 
and  then  defines  two  weighting  matrices  in  preparation  for  the  study  of  a  linear  quadratic  regulator 
problem.  Rather  than  inputting  the  matrices  from  the  keyboard,  the  RDF  (Read  Data  File) 
operator  is  used.  Since  the  matrix  E  is  equal  to  the  identity  for  this  problem,  the  terminal 
constraint  (2.1.14)  reduces  to  K  =  Q  which  is  performed  in  step  13  below. 

1  :Linear  Quadratic  Regulator  Problem 

2  ; 

3  ;System  and  Weighting  Matrix  Identification 

4  :  —  System  Definition  — 

5  (rdf)=a,b,c 

|  6  a.b.c(out)* 

7  a.b.c(syst)= 

8  :  —  Weighting  Matrices  Definition  — 

9  (rdf)=q.r 

10  q,r(out)= 

11  q.r(lq)= 

1  12  ^Initialization 

13  q(mcp)=k 

14  l(dsc)=one 

15  ; 

16  "Enter  the  total  number  of  iterations  to  perform.  (A  scalar)" 

17  (inp)=num 

|  18  ; 

19  ;  Main  Loop 

'  20  Z:k(out)= 

21  k(dre)-knew 

|  22  knew(mcp)-k 

23  num.one(-)«num 

24  num(if)=Z 

!  25  : 

1  26  ;  Done.  Print  out  final  k. 

27  k(out)- 

|  28  (stop)* 

|  If  this  program  is  run  for  five  iterations  using  the  second-order  system  of  Section  4.1.  the  following 

|  output  results. 


* 


/• 


i 

I 


.V 
»  _ 


>  ;Linear  Quadratic  Regulator  Problem 

>  ; 

>  iSvstem  and  Weighting  Matrix  identification 


V 


> :  —  System  Definition 
>(rdf)**a,b.c 


Enter  name  of  the  Data  File  (DF)  for  matrix  a  >syst 

Opening  file  named  :  syst.DF 

Reading  array  named  :  a 

Reading  array  named  :  b 

Reading  array  named  :  c 

>a,b,c(out)« 

a 

1.000  2.000 
3.000  4.000 

b 

1.000 

0. 

c 

0.  1.000 
>a.b.c(syst)= 

> ;  —  Weighting  Matrices  Definition  — 

>(rdf)=q.r 

Enter  name  of  the  Data  File  (DF)  for  matrix  q  >  lq 
Opening  file  named  :  lq.DF 
Reading  array  named  :  q 
Reading  array  named  :  r 

>q.r(out)“ 

q 

1.000  0. 

0.  2.000 

r 

1.000 

>q.r(lq)= 

>  initialization 
>q(mcp)=k 

>  l(dst  )=one 
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>  "Enter  the  total  number  of  iterations  to  perform.  (A  scalar)" 
Enter  the  total  number  of  iterations  to  perform.  (A  scalar) 

>(inp)-=num 
***  Matrix  num  *** 

Enter  the  dimensions  of  this  matrix.  >1.1 
Enter  the  scalar  :  num  >  5 
> ; 

> ;  Main  Loop 
>Z:k(out)- 
k 

1.000  0. 

0.  2.000 

>k(dre)=knew 

>knew(mcp)=k 

>  num.one(-)=num 
>num(if)=Z 
>Z:k(oul)= 

k 

19.500  25.000 
25.000  36.000 

>  k(dre)=knew 

>  knew(mcp)=k 
>num.one(-)=num 

>  num(if  )=Z 
>Z:k(out)= 

k 

58.878  80.244 
80.244  113.512 

>  k(dre)=knew 

>  knew(  mcp)=k 


>  num.one(-)«=num 
>num(if)-Z 
>Z:k(out)= 

k 

63.804  87.075 
87.075  122.984 

>k(dre)=knew 

>knew(mcp)«k 

>  num.one(-)=num 
>num(if)=Z 
>Z:k(out)= 

k 

63.918  87.234 
87.234  123.208 

>  k(dre)=knew 

>  knew(mcp)=k 
>num.one(-)=*num 
>num(if)=Z 

> ; 

>  ;  Done.  Print  out  final  k. 

>  k(out)= 

k 

63.922  87.240 
87.240  123.216 

>(stop)= 


Notice  the  speed  with  which  the  Riccati  iterations  converge.  The  value  of  k  is  settled  to  three 
significant  digits  in  only  five  iterations. 
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4.4  L— A— S  Operator  GAME 


•  < 


The  L-A-S  operator  named  GAME  is  used  to  define  a  linear  shift-invariant  descriptor  game. 
Often,  it  is  the  first  operator  issued  to  the  L-A-S  interpreter  when  a  linear  quadratic  descriptor 
Nash  game  problem  is  being  studied. 

L  —A  —S  Help  —File  Description 

GAME  -  2-player  GAME  description 
Syntax  :  A.  Bl.  B2.  Cl.  C2  [,E  ]  (GAME)  = 

Input  Data  :  A  [N.N]  .  Bl  [N.Pl]  .  B2  [N.P2]  . 

Cl  [M1.N]  .  C2  [M2.N]  .  E  [N.N] 

Matrix  E  is  optional. 

Options  :  E,  L.  T 

Description  :  Identifies  the  following  discrete-time  game  : 


£ 


© 


5 

t* 


3 


E  x(k+l)  =  A  x(k)  +  BjUjCk)  +  B2u 2(k) 
y  !(k)  =  C[Jc(k) 
y  2(k)  =  C2  x  (k) 


y 

A' 


Note  :  If  E  is  omitted,  it  is  assumed  to  be  the  identity  matrix. 


Example  of  Usage 

An  example  of  the  usage  of  operator  GAME  is  given  below.  The  following  L-A-S  program 

segment  identifies  a  discrete-time  descriptor  game. 

;  Descriptor  Game  Identification 
(rdf)=a.bl.b2,cl.c2 
a.bl.b2.cl.c2(out)= 
a.bl  .b2.c  1  .c2(game)= 

If  this  program  is  executed  for  a  third-order  game  problem,  the  following  output  would  result. 


>  :  Descriptor  Game  Identification 

>  (rdf  )=a.bl  ,b2.cl  ,c2 

Enter  name  of  the  Data  File  (DF)  for  matrix  a  >ab 

Opening  file  named  :  ab.DF  I 

Reading  array  named  :  a 
Reading  array  named  :  bl 

Reading  array  named  :  b2  .'•! 
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fit  »»*. 


E 


Enter  name  of  the  Data  File  (DF)  for  matrix  cl  >cc 
Opening  file  named  :  cc.DF 
Reading  array  named  :  cl 
Reading  array  named  :  c2 

>  a.bl  ,b2.cl  ,c2(out)= 


0.435  -1.401  -0.896 
-0.172  -0.569  1.391 
-1.655  0.008  0.134 

bl 

1.000  2.000 
3.000  4.000 
5.000  6.000 


i 


1.000 

'm 

1.000 

0. 

* 

cl 

1.000  0. 

0. 

0.  1.000 

0. 

| 

0.  0. 

1.000 

c2 

v: 

1.000  0. 

0. 

0.  1.000  0. 

0.  0.  1.000 

>  a  .b  1  .b2  ,c  1  .c2(  game)= 

4.5  L— A— S  Operator  LQNG 

The  L-A-S  operator  named  LQNG  is  used  to  iterate  two  coupled  discrete-time  Riccati 
equations. 

L—A—S  Help— File  Description 

LQNG  -  Linear  Quadratic  Nash  Game  (Coupled  Riccati  Iterations) 


Syntax 
Input  Data 
Output  Data 

Options 

Description 


K1 .  K2  (LQNG)  =  KIN.  K2N  [.  I  I  [.  F2]  ] 

K 1  [N.N]  .  K2  [N.N] 

KIN  [N.N]  .  K2N  [N.N]  .  FI  [Pl.N]  .  F2  [P2..N] 
Arravs  FI  and  F2  are  optional. 

E.  L.  T 

Assuming  that  operators  GAME  and  LQ  have  been  issued 
this  operator  performs  one  iteration  of  two  coupled 
discrete-time  Riccati  equations.  KIN  and  K2N  are  the 


\  ,V.\  • 


new  Riccati  gain  matrices.  FI  and  F2  (if  provided) 
are  the  corresponding  feedback  matrices. 


Example  of  Usage 

An  example  of  the  usage  of  operator  LQNG  is  given  below.  Consider  the  discrete-time  game 
defined  in  Section  4.4.  The  following  segment  is  an  L-A-S  program  that  identifies  the  same  system 
and  then  defines  two  weighting  matrices  in  preparation  for  the  study  of  a  linear  quadratic  Nash 
game.  Rather  than  inputting  the  matrices  from  the  keyboard,  the  RDF  (Read  Data  File)  operator  is 
used.  Since  the  matrix  E  is  equal  to  the  identity  for  this  problem,  the  terminal  constraint  (2.1.14) 
reduces  to  K x  -  (2,  which  is  performed  in  step  10  below. 

1  :  Linear  Quadratic  Nash  Game 

2  (rdf)=a.bl.b2.cl.c2 

3  a.bl.b2.cl.c2(out)= 

4  a.bl.b2.cl.c2(game)= 

5  (:rdf)=rl.r2.sl.s2 

6  rl.r2,sl.s2(out)= 

7  cl(t).sl(*).cl(*)=ql 

8  c2(t).s2(*).c2(*)=q2 

9  ql.rl.q2.r2(lq)= 

10  ql.q2(mcp)=kl.k2 

11  ; 

12  l(dsc)=one 

13  "Enter  the  total  number  of  stages  in  this  game." 

14  (inp)=ii 

15  :  Main  Loop 

16  a:kl.k2(out)= 

17  kl.k2(lqng)=kln,k2n 

18  kln.k2n(out)= 

19  kln.k2n(mcp)»<kl,k2 

20  ii,one(-)=ii 

21  ii(  if  )=a 

If  this  program  is  run  for  five  iterations  using  the  third-order  system  of  Section  4.4.  the  following 
output  results. 


>  ;  Linear  Quadratic  Nash  Game 

>  ( rdf  )=a,bl  .b2.cl  ,c2 

Enter  name  of  the  Data  File  (DF)  for  matrix  a  >ab 
Opening  file  named  :  ab.DF 
Reading  array  named  :  a 
Reading  array  named  :  hi 


77 


HH 


Reading  array  named  :  b2 

Enter  name  of  the  Data  File  (DF)  for  matrix  cl 
Opening  file  named  :  cc.DF 
Reading  array  named  :  cl 
Reading  array  named  :  c2 

>  a.bl  .b2  ,cl  .c2(out)= 


0.435  -1.401  -0.896 
-0.172  -0.569  1.391 
-1.655  0.008  0.134 

bl 

1.000  2.000 
3.000  4.000 
5.000  6.000 

b2 

1.000 

1.000 

0. 

cl 

1.000  0.  0. 

0.  1.000  0. 

0.  0.  1.000 

c2 

1.000  0.  0. 

0.  1.000  0. 

0.  0.  1.000 

>  a.bl  ,b2.cl  .c2(game)= 

>  ( rdf  )=r  1  .r2  .s  1  .s2 

Enter  name  of  the  Data  File  (DF)  for  matrix  rl  >rs 

Opening  file  named  :  rs.DF 

Reading  array  named  :  rl 

Reading  array  named  :  r2 

Reading  array  named  :  si 

Reading  array  named  :  s2 

>  rl  ,r2.sl  ,s2(out)= 

rl 

1.000  0. 

0.  1.000 


r2 

1 .000 
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si 

1.000  0.  0. 

0.  1.000  0. 

0.  0.  1.000 

s2 

1.000  0.  0. 

0.  1.000  0. 

0.  0.  1.000 

>cl(t),sl(*),cl(*)=ql 

#l.sl(*).cl(*)~ql 

#2.cl(*)«ql 

>c2(t).s2(*).c2(*)-q2 
#1  (s2(*),c2(*)«q2 
#2.c2(*)»q2 

>  q  1  ,r  1  .q2  ,r2(  lq  )= 

>  q  1  ,q2(  mcp)=k  1  .k2 

>: 

>  l(dsc)=one 

>  "Enter  the  total  number  of  stages  in  this  game." 
Enter  the  total  number  of  stages  in  this  game. 

>(inp)=ii 
***  Matrix  ii  *** 

Enter  the  dimensions  of  this  matrix.  >1.1 
Enter  the  scalar  :  ii  >5 


> :  Main  Loop 


>a:kl  .k2(out)= 
kl 

1.000  0.  0. 

0.  1.000  0. 

0.  0.  1.000 

k2 

1.000  0.  0. 

0.  1.000  0. 

0.  0.  1.000 


>kl.k2(lqng)=kln.k2n 


i 
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>kln,k2n(out)= 


1 

kin 

w 

1.236  -0.229 

0.003 

-0.229  1.569 

0.817 

£ 

£ 

0.003  0.817 

3.291 

k2n 

w 

1.393  -0.432 

0.090 

-0.432  1.771 

0.709 

0.090  0.709 

3.235 

<'« 

& 

>  kln.k2n(mcp)=kl  ,k2 

%*’ 

>ii.one(-)=ii 

>/ 

*  *’ 

>ii(if)=a 

>a:kl.k2(out)= 

■ 

kl 

v 

1.236  -0.229  0.003 

$ 

-0.229  1.569  0.817 

w 

0.003  0.817  3.291 

1 

k2 

a 

1.393  -0.432  0.090 
-0.432  1.771  0.709 
0.090  0.709  3.235 

f i 

>kl  ,k2(lqng)=kln.k2n 

% 

>kln.k2n(out)= 

kin 

>■ 

1.242  -0.238  -0.009 

-0.238  1.715  1.192 

-0.009  1.192  4.282 

k2n 

1.402  -0.418  0.159 

-0.418  1.971  1.302 

£ 

0.159  1.302  5.081 

■  ■ 

>kln,k2n(mcp)=kl.k2 

> 

«'*■ 

>  ii.one(-)=ii 

| 

>  ii(  if  )=a 

>a:kl  ,k2(out)= 

A 

kl 
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1.242  -0.238  -0.009 
-0.238  1.715  1.192 
-0.009  1.192  4.282 


k2 

1.402  -0.418  0.159 
-0.418  1.971  1.302 
0.159  1.302  5.081 

>  k  1  ,k2(  lqng)«k  1  n  ,k2n 


>kln.k2n(oui)= 


kin 

1.247  -0.238  0.007 
-0.238  1.710  1.181 
0.007  1.181  4.304 

k2n 

1.403  -0.416  0.168 
-0.416  1.977  1.321 
0.168  1.321  5.155 

>  k  1  n  ,k2n(mcp)-=k  1  .k2 

>  ii.one(-)=ii 


>ii(if)=a 


>a:kl,k2(out)= 

kl 

1.247  -0.238  0.007 
-0.238  1.710  1.181 
0.007  1.181  4.304 

k2 

1.403  -0.416  0.168 
-0.416  1.977  1.321 
0.168  1.321  5.155 

>  kl  ,k2(lqng)=k  ln.k2n 

>kln,k2n(out)= 


kin 

1.247  -0.238  0.006 
-0.238  1.714  1.192 
0.006  1.192  4.330 

k2n 

1.403  -0.416  0.169 
-0.416  1.976  1.320 


81 


0.169  1.320  5.155 

>  k  1  n  ,k2n(  mcp  )=k  1  .k2 

>  ii.one(-)=ii 

>  ii(if  )— a 
>a:kl,k2(out)» 

kl 

1.247  -0.238  0.006 
-0.238  1.714  1.192 
0.006  1.192  4.330 

k2 

1.403  -0.416  0.169 
-0.416  1.976  1.320 
0.169  1.320  5.155 

>  k  1  .k2(  lqng)=k  1  n.k2n 
>kln.k2n(out)= 

kin 

1.247  -0.238  0.006 
-0.238  1.714  1.192 
0.006  1.192  4.330 

k2n 

1.403  -0.416  0.169 
-0.416  1.976  1.320 
0.169  1.320  5.156 

>  k  1  n,k2n(  mcp)=k  1  .k2 
>ii.one(-)=ii 

>  ii(  if  )=a 


4.6  L  -A— S  Operator  MLTR 

The  L-A-S  operator  named  MLTR  is  used  to  iterate  two  coupled  discrete-time  Riccati 
equations  where  multirates  are  involved. 


L—A—S  Help— File  Description 

MLTR  -  MuLTiRate  nash  game  (Multirate  Coupled  Riccati  Iterations) 
Syntax  :  K1.K2.  K.  X  (MLTR)  -  KIN.  K2N  [.  FI  [.  F2]  ] 
Input  Data  :  Kl  [\.\]  .  K2  [N.N]  .  K  [l  .1  ]  .  N[l.l] 
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Output  Data:  KIN  [N.N]  .  K2N  [N.N]  .  F1[P1,N]  .  F2  [P2.N] 

Arrays  FI  and  F2  are  optional. 

Options  :  E,  L.  T 

Description  :  Assuming  that  operators  GAME  and  LQ  have  been  issued. 

this  operator  performs  one  iteration  of  two  coupled 
discrete-time  Riccati  equations  where  multirates 
are  involved.  K  is  the  current  time  instant  and 
N  is  the  multirate  parameter.  KIN  and  K2N  are  the 
new  Riccati  gain  matrices.  FI  and  F2  (if  provided) 
are  the  corresponding  feedback  matrices. 


Example  of  Usage 

An  example  of  the  usage  of  operator  MLTR  is  given  below.  Consider  a  third-order  discrete¬ 
time  system  that  is  being  controlled  by  two  computers  operating  at  different  speeds.  The 
computers  are  decentralized  in  that  they  are  located  at  different  physical  places.  They  control  the 
plant  through  a  telephone  hookup.  DM1  is  equipped  with  a  1200  baud  modem,  but  DM2  has  only 
a  300  baud  modem.  Thus.  N  —  4  for  this  problem.  Because  of  the  interface  between  the  plant  and 
the  controllers,  the  input  of  DM2  is  the  all-digital  control  policy  (2.2.3).  The  following  L-A-S 
program  is  used  to  study  this  multirate  game. 

1  :  Linear  Quadratic  Multirate  Nash  Game 

2  (rdf)=a.bl.b2.cl.c2 

3  a.bl.b2,cl,c2(out)= 

4  a.bl.b2.cl.c2(game)= 

5  (rdf)=sl.rl.s2.r2 

6  cl(t).sl(*).cl(*)=»ql 

7  c2(t).s2(*).c2(*)=q2 

8  ql.rl .q2.r2(out)= 

9  ql.rl.q2,r2(lq)= 

10  ; 

11  ql.q2(mcp)=kl.k2 

12  l(dsc)=one 

13  Enter  the  total  number  of  stages  in  this  game" 

14  (inp)=ii 

15  "Enter  the  multirate  parameter.  N" 

16  (inp)=N 

1 7  ;  Main  Loop 

18  a:kl.k2(out)= 

19  kl  ,k2.ii.N(mltr)=kln.k2n.fl  ,f2 

20  1 1  ,f2(out.e)= 

21  kln.k2n(mcp)=kl  ,k2 

22  ii.one(-)=ii 

23  ii(if )=a 

24  ; 

25  kl  .k2.ii.N(mltr)=k  ln.k2n.fl  .12 


o 


26  fl,f2(out.e)« 


If  this  program  is  run  for  six  iterations  using  an  arbitrary  third-order  system,  the  following  output 


results. 


> ;  Linear  Quadratic  Multirate  Nash  Game 
Xrdf)=a,bl.b2.cl.c2 

Enter  name  of  the  Data  File  (DF)  for  matrix  a  >  ACCl 

Opening  file  named  :  ACCl.DF 

Reading  array  named  :  a 

Reading  array  named  :  bl 

Reading  array  named  :  b2 

Reading  array  named  :  cl 

Reading  array  named  :  c2 

>a.bl  ,b2.cl,c2(out)= 


0.435  -1.401  -0.896 
-0.172  -0.569  1.391 
-1.655  0.008  0.134 

bl 

1.000 

0. 

1.000 

b2 

1.000 

1.000 

0. 

cl 

1.000  0.  0. 

0.  1.000  0. 

0.  0.  1.000 

c2 

1.000  0.  0. 

0.  1 .000  0. 

o.  0.  1.000 

>  a.bl  ,b2  ,c  1  .c2(game)= 

>  (rdf  )=sl  ,rl  ,s2.r2 

Enter  name  of  the  Data  File  (DF)  for  matrix  si  >  ACC2 
Opening  fi’e  named  :  ACC2.DF 
Reading  array  named  :  si 


o' ^ -o'- .  ■ 


a 


Reading  array  named  :  rl 
Reading  array  named  :  s2 
Reading  array  named  :  r2 

>cl(t).sl(*).cl(*)=ql 

#l.sl(*).cl(*)-ql 

#2.cl(*)-ql 

>  c2(t),s2(*).c2(*)=q2 
#l.s2(*).c2(*)-q2 
#2.c2(*)=q2 

>  q  1  .rl  ,q2  ,r2(out  )= 

qi 

1.000  0.  0. 

0.  1.000  0. 

0.  0.  1.000 

rl 

1.000 

q2 

1.000  0.  0. 

0.  1.000  0. 

0.  0.  1.000 

r2 

1.000 


>ql.rl.q2.r2(lq)= 


>ql  .q2(mcp)=kl  .k2 

>  l(dsc)=one 

>  "Enter  the  total  number  of  stages  in  this  game" 
Enter  the  total  number  of  stages  in  this  game 

>  ( in  p)=ii 

***  Matrix  ii  *** 

Enter  the  dimensions-  of  this  matrix.  >1.1 
Enter  the  scalar  :  ii  >6 

>  "Enter  the  multirate  parameter,  \" 

Enter  the  multirate  parameter,  N 


*  *  *.•'*»'*  V  •'*  *  " 


>v*  .>\ v .v. 
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>(inp)-N 
***  Matrix  N  *** 

Enter  the  dimensions  of  this  matrix.  >1.1 
Enter  the  scalar  :  N  >4 
> ;  Main  Loop 
>a:kl.k2(out)- 
kl 

1.000  0.  0. 

0.  1.000  0. 

0.  0.  1.000 

k2 

1.000  0.  0. 

0.  1.000  0. 

0.  0.  1.000 

>  kl  .k2.ii,N(mltr)=k  ln.k2n.f  1  ,f2 
>fl.f2(out.e)= 

fl 

-4.06588e-01  -4.64582e-01  -2.541 76e-01 
f2 

0.00000e+00  0.00000e+00  0.00000e+00 

>  k  1  n  .k2n(  mcp)=k  1  .k2 

>  ii.one(-)=ii 
>ii(if)=a 
>a:kl.k2(out)= 

kl 

3.461  -1.091  -1.160 
-1.091  2.640  0.112 
-1.160  0.112  3.562 


k2 

3.295  -1.280  -1.264 
-1.280  2.424  -0.006 


-4.91982e-01  -4.64462e-01  -5.44093e-01 


f2 

0.00000e+00  O.OOOOOe+OO  0.00000e+00 
>kln,k2n(mcp)=kl.k2 
>ii.one(-)=ii 
>ii(if)=a 
>a:kl.k2(out)= 
kl 

12.002  -5.782  -7.172 
-5.782  5.706  2.611 
-7.172  2.611  10.304 

k2 

11.760  -6.010  -7.441 
-6.010  5.214  2.712 
-7.441  2.712  9.552 

>kl  ,k2.ii,N(mltr)=kln.k2n.f  l.f2 

>f  l.f2(out.e)= 

fl 

-6.281 13e-01  -3.24123e-01  -7.37706e-01 
f2 

1.86462e+00  -1.22783e+00  -1.03105e+00 

>  kln.k2n(mcp)=kl.k2 

>  ii.one(-)=ii 

>  ii(  if  )=a 


>a:kl.k2(out)= 

kl 

23.862  -8.100-23.483 
-8.100  4.268  9.051 
-23.483  9.051  27.654 

k2 

23.198  -9.483  -21.761 
-9.483  5.342  8.857 
-21.761  8.857  23.593 


>  k  1  .k2.ii..\’(mltr)=k ln,k2n,f  1  ,f2 
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■  >  -s. 


>fl,f2(out.e)» 


-1.24307e+00  -1.87479e-01  2.77446e-01 


0.00000e+00  0.00000e+00  0.00000e+00 


>  kln.k2n(mcp)-kl  ,k2 
>ii.one(-)-ii 

>  ii(  if  )=— a 
>a:kl.k2(out)= 


112.932  -61.747  -77.912 
-61.747  36.568  42.550 
-77.912  42.550  57.697 


107.254  -57.781  -79.740 
-57.781  33.065  42.170 
-79.740  42.170  63.870 

>kl.k2.ii.N(mltr)=kln,k2n.fl.f2 

>fl.f2(out.e)=*= 


3.28894e+00  -2.42437e+00  -3.84678e+00 


0.00000e+00  0.00000e+00  0.00000e+00 


>  k  1  n.k2n(mcp)=k  1  ,k2 

>  ii.one(-)=ii 

>  ii(  if  )=a 
>a:kl.k2(out)= 


1.56022e+02  -8.05226e+01  -1.31445e+02 
-8.05226e+01  4.46216e+01  6.68516e+01 

-1.31445e+02  6.68516e+01  1.18129e+02 


2.014()le+02  -1.03781e+02  -1.59092e+02 
-1.03781e+02  5.55282e+01  8.09325e+01 

-1.59092e+02  8.09325e+01  1.30365e+02 


S 

s 

.» 

.4 
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>kl.k2.ii.N(mltr)=kln.k2n.fl.f2 
>f  l,f2(out,e)*= 
fl 

2.86026e+00  -2.18331e+00  -3.49269e+00 
f2 

0.00000e+00  0.00000e+00  0.00000e+00 

>  k  1  n  ,k2n(  mcp)=k  1  ,k2 

>ii.one(-)=ii 

>ii(if)=a 

> ; 

>kl,k2,ii.N(mltr)=kln.k2n,fl.f2 

>fl.f2(out.e)= 

fl 

-4.35703e-01  -4.13821e-01  -9.23808e-01 
f2 

3.93805e+00  -2.10267e+00  -3.12784e+00 
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CHAPTER  5 
CONCLUSIONS 

This  dissertation  studies  the  computational  aspects  of  iterating  two  coupled  discrete-time 
Riccati  equations.  These  equations  arise  as  a  result  of  solving  LQ  descriptor  Nash  games.  The 
presence  of  coupling  in  the  Riccati  equations  complicates  the  iteration  process.  This  work  devises  a 
method  which  removes  the  coupling  in  a  numerically  robust  manner.  Then  algorithms  are 
engineered  that  compute  the  quantities  needed  to  iterate  the  Riccati  equations.  Every  opportunity 
is  taken  to  exploit  the  properties  of  matrices  (e.g.,  positive-definiteness)  that  enter  into  the 
calculations  so  as  to  obtain  a  savings  in  computation.  The  algorithms  are  coded  and  the  coupled 
Riccati  software  is  integrated  into  the  L-A-S  CACSD  language.  The  novelty  of  this  work  is  two¬ 
fold.  First,  a  new  problem  is  formulated,  solved  and  the  solution  procedure  is  implemented  as 
computer  code.  Second,  numerical  theory  and  software  are  combined  under  the  heading  of  CACSD 
to  yield  a  package  that  allows  others  10  solve  single  and  coupled  Riccati  equations. 

The  numerical  issues  associated  with  iterating  coupled  discrete-time  Riccati  equations  are  the 
key  focus  of  this  thesis.  The  software  developed  for  the  iteration  task  is  coded  in  FORTRAN  and 
makes  extensive  use  of  the  LINPACK  library.  A  structured  programming  approach  has  produced 
low-level  algorithms  that  are  very  modular  and  extremely  efficient.  The  final  result  is  a  set  of  six 
new  L-A-S  operators  that  are  collectively  capable  of  iterating  single  and  coupled  discrete-time 
Riccati  equations. 

Although  the  main  contribution  of  this  thesis  is  the  software  engineering  of  the  coupled 
Riccati  problem,  there  are  several  theoretical  advancements  which  add  breadth  to  the  work.  The 
most  important  of  these  are  the  existence  and  convergence  theorems  which  define  the  iteration 
behavior  for  both  finite-horizon  and  infinite-horizon  problems.  Of  less  relevance  but  not 
significance  is  the  development  of  descriptor-variable  dynamic  games.  In  addition  to  theoretical 
extensions,  there  are  physical  and  numerical  advantages  associated  with  descriptor  game 
formulations.  Multirate  I.Q  Nash  games  were  heretofore  unposed. 
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There  are  several  directions  for  future  research.  For  example,  the  assumption  that  E  is 
nonsingular  could  be  relaxed.  Then  uniqueness  of  the  Riccati  iterates  is  lost.  Less  restrictive 
conditions  insuring  convergence  for  infinite- horizon  problems  could  be  sought.  A  contraction 
mapping  argument  could  be  attempted  for  the  case  of  unstable  plants.  Lyapunov-type  stability 
results  could  be  applied  to  the  coupled  Riccati  maps.  A  completely  different  research  topic  would 
be  the  determination  of  Leader-Follower  strategies  for  descriptor  games. 
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APPENDIX  A 


SOFTWARE  LISTINGS 


This  appendix  contains  the  FORTRAN  computer  codes  used  for  a  selected  number  of  the  low- 
level.  coupled  Riccati  algorithms.  In  alphabetical  order,  the  algorithms  are  MLTPLY.  PSICOM  and 
QDFORM. 


Algorithm  MLTPLY 


************************* 

**  M  L  T  P  L  Y  ** 

************************* 

*  * 

**  Given  the  double  precision  matrices  A  and  B,  this  subroutine 
**  computes  and  returns  the  variable  C  defined  as  the  product 
*  * 

**  C  =  A  *  B  ,  if  the  logical  variable  ATFLAG  is  False,  or 

*  * 

*  *  *p 

**  C  »  A  *  B.  if  the  logical  variable  ATFLAG  is  True. 

*  * 

**  where  A  is  an  NxM  arbitrary  matrix  (for  ATFLAG  »  False). 

**  is  an  MxN  arbitrary  matrix  (for  ATFLAG  -  True). 

*  *  B  i s  a  n  Mx  P  arbitrary  ma  t  r  i  x  , 

**  and  C  is  the  resulting  NxP  matrix. 


SUBROUTINE  MLTPLY  ( A . B . C .N .M.  P . ATFLAG ) 


DOUBLE  PRECI SION 

INTEGER 

LOGICAL 

► 

INTEGER 


A  (  1  ) 

N  .  M 
ATFLAG 


B  (  1  ) 
P 


C  (1) 


I.I1.I2.I3.J.K 


'  Branch  depending  on  the  value  of  ATFLAG. 

t 

IF  (ATFLAG)  GO  TO  2 

1  Compute  A  *  B  and  store  in  C 

t 

DO  1  J  =  1 . P 
13  =  ( J  -  1  )  *N 

DO  1  1 =1 ,N 

11  =  I  -  N 

12  =  ( J -  1  )  *M 

13  =13+1 
C  (13)  =  0.0 


V  V  v.v.  .■ 


im 


m 


m 


TwYV  "in” 


vi?  «•.-  *.j<  «j.  ■-'.'-r-T^v'-.1-"-  ,i.v  ■•  '\ 
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C  ** 

DO  1  K-l .M 

11  -  1 1  +  N 

12  -12+1 

C  (13)  -  C  (13)  +  A  (II)  *  B  (12) 

1  CONTINUE 
RETURN 

C  ** 

C  *  *  T 

C  **  Compute  A  *  B  and  store  in  C. 

C  ** 

2  CONTINUE 

DO  3  J-l .P 
13  -  ( J-l  )*N 

DO  3  1=1 .N 

11  -  (  I  -  1  )  *M 

12  -  ( J -  1  )  *M 

13  =13+1 
C  (13)  -  0.0 

c  ** 

DO  3  K=1 .M 

11  =11+1 

12  =12+1 

C  (13)  =  C  (13)  +  A  (II)  *  B  (12) 

3  CONTINUE 
RETURN 
END 


$ 


UUUUUUUUUUUUUUUUUUUUUUUUU'J  uooooo 
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Algorithm  PS1COM 


*********************** 

P  S  1  C  O  M  ** 
***************** 


Given  B.  KK  and  R.  this  subroutine  computes  and  returns 
the  variable  named  PS  I  defined  as 

-1  T 

PS  I  -  GAMMA  *  B  *  KK 

T 

where  GAMMA  -  R  +  B  ♦  KK  *  B 
B  i s  a  NxP  matrix. 

KK  is  a  NxN .  po s i t i ve - s emi  d e f i n  i  t e  matrix, 
and  R  is  a  PxP.  positive  matrix. 

Note  It  is  assumed  that  P  is  less  than  or  equal  to  N! 

Several  other  quantities  are  calculated  and  returned  for 
possible  later  use.  These  quantities  include 
T 

W  -  B  *  KK  .  GAMMA  .  XK .  and  XG  where  XG  is  the  Cholesky 

T 

factor  of  GAMMA .  That  is.  GAIVMA  ■  XG  *  XG . 

Likewise  XK  is  the  Cholesky  factor  of  KK . 

SUBROUTINE  PSICCM  ( B . GAMMA . KK , R  W.  XG . XK .N . P . PS  I ) 


GAMMA  contains  enough  ro  om  for  1  double  precision  PxP  ma  t  r  i  x  . 
W  contains  enough  room  for  1  double  precision  NxN  matrix. 

XG  contains  enough  room  for  1  double  precision  PxP  matrix. 

XK  contains  enough  room  for  1  double  precision  NxN  matrix. 


C 

C 

c 

c 

c 

c 

c 

c 

c 

c 


DOUBLE  PRECISION 
DOUBLE  PRECISION 
INTEGER 

INTEGER 

DATA 


B  (  1  )  GAMMA  (  1  ) 

W  ( 1  )  .  XG  ( 1 )  . 

N  .  P 

l  .  I  1  .  J  .  K  .  L 

USER  76/ 


.  KK  (  1  )  .  R  ( 1 ) 

XK  (  1  )  .  PS  I  (  1  ) 

USER 


*  *  J 

*  *  Compu t  e  B  *  KK  and  store  in  W  using  algorit  hm  MLTPLY . 

*  * 


CALL  MLTPLY  (  B  .  KK  ,W.  P  .  N  .  N  .  TRUE  .  ) 

*  * 

*  *  J 

**  Compute  GAMMA  =  R  +  B  *  KK  *  B  using  algorithm  QDFORM. 
**  Note  Since  XG  is  only  computed  later,  it  is  used  here  as  a 


Cl  O'  Ci  uiOOOOO 


PIIWWJPP..I.LJ  J  J  J 


7.T!p.jLij|i|  (j  "1  T*  -> -J>  v  v  .* 


C 

c 

C 

C 

c 

c 

c 


3 

c 


4 


7 
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**  dummy  variable  for  temporary  storage  by  QDFORM. 

**  Also.  CHOFLG  -  False.  CFLAG  -  True,  and  ADDFLG  -  True. 

*  * 

CALL  QDFORM  (KK.B.R, GAMMA. XK.XG.N.P.  .FALSE.  .  .TRUE.  .  TRUE.) 

*  * 

» *  7 

**  Factor  GA\MA  =  XG  *  XG  wh ere  XG  is  upper  triangular. 

*  * 

K  -  P  *  P 

DO  3  1-1 .K 

XG  (  I  )  -  GAMV1A  (  I  ) 

CONTINUE 

*  * 

CALL  DPOFA  (XG.P.P.I) 

IF  ( I  . EQ .  0 )  GO  TO  5 

WRITE  (USER. 4) 

FORMAT  (/•  ERROR  :  PSICGM  -  GAIVMA  is  not  positive  definite!’) 
RETURN 

«  « 

**  Compute  PSI  by  solving  the  set  of  linear  equations 

*  *  7 

**  GAMMA  *  PSI  -  B  *  K  =  W. 

*  * 

CONTINUE 
II  =  -P 

DO  7  J  =  1 ,N 

**  Copy  Jth  column  of  W  to  Jth  column  of  PSI. 

II  -  II  +  P 

L  =11 

DO  6  K= 1 . P 
L  =  L  +  1 

PSI  (L)  =  W  (L) 

CONT I NUE 

**  Compute  Jth  column  of  PSI. 

CALL  DPOSL  ( XG .P.P, PS  1(11  +  1)) 

CONT I NUE 
RETURN 
END 


>: 

y. 


?: 

& 


r. 


i 
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Algorithm  QDFORM 


*********************** 


QDFORM 


Given  A.  B  and  optionally  C.  this  subroutine  computes  and 
returns  one  of  the  following  quadratic  forms 


B  *  A  *  B 


C  +  B  *  A  *  B 


D  -  C  -  B  *  A  *  B 

whe  r  e  A  is  a  NxN .  (positive -definite)  ma  t  .  *  x  , 

B  is  a  NxP  arbitrary  matrix, 
and  C  is  a  PxP  symmetric  matrix. 

It  is  assumed  that  P  is  less  than  or  equal  to  N! 

SUBROUTINE  QDFORM  (A.B.C.D.X.Y.N.P. CHOFLG . CFLAG . ADDFLG ) 

Note  that  the  arrays  X  and  Y  must  be  provided  by  the  calling 
routine  as  temporary  storage  for  the  calculation.  Each 
must  contain  enough  room  for  1  double  precision  NxN  matrix. 

This  routine  has  three  flags  which  govern  the  calculation  of  a 
quadratic  fo  rm . 

If  CHOFLG  is  True,  then  the  Cholesky  factor  of  A  is  already 
available  and  has  been  passed  in  X.  Consequently,  the  call 
to  L INPACK  routine  DPOFA  is  skipped.  If  CHOFLG  is  False, 
then  the  Cholesky  factor  of  A  is  assumed  to  be  unavailable. 

If  CFLAG  is  True,  then  the  matrix  C  is  to  be  included  in  the 
initialization  statement  of  the  matrix  multiply  DO-Loop. 

If  CFLAG  is  False,  then  C  is  never  referenced. 

Assuming  CFLAG  -  True,  then  ADDFLG  is  consulted  to  determine 
if  an  addition  or  subtraction  is  to  be  performed  in  the 
matrix  multiply  DO-Loop.  .ADDFLG  =  True  means 
T 

C  +  B  *  A  *  B  is  computed  ADDFLG  =  False  means 

T 

C  -  B  *  A  *  B  is  computed.  If  CFLAG  =  False,  then 

.ADDFLG  is  never  referenced. 


DOUBLE  PRECISION 


A  (  1  ) 


B  (  1  )  .  C  (  1  ) 


D  (  1  ) 


Atm 


A.'uV.^VoVvVvS'>-jV>,v 


2a 


c 

c 


c 

c 

c 


c 

c 

c 

c 


1 

c 


c 

c 

c 

c 

3 


4 

c 

c 

c 


c 

c 

c 

c 

c 
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DOUBLE  PRECISION 

INTEGER 

LOGICAL 


X  (1  ) 

N  .  P 
CHOFLG 


INTEGER 


DATA 


I  .  I  1  . 
USER  /6 / 


Y  (1) 
CFLAG 
12.13 


ADDFLG 
14  .  J  .  K 


USER 


**  Skip  the  Cholesky  factorization  of  A  if  already  available. 

*  * 


IF 


(CHOFLG)  GO  TO  3 


*  * 
*  * 
*  * 
*  * 


Factor  A 


T 

X  *  X 


wh ere  X  is  upper  triangular. 


*  * 
*  * 


K  -  N  *  N 

DO  1  1  =  1  .K 

X  (  I  )  =  A  (  I  ) 

CONTINUE 

s 

CALL  DPOFA  ( X . 

N.N.  I  ) 

IF  ( I  . EQ .  0 ) 

WRITE  (USER. 2) 

GO  TO  3 

FORMAT  ( / •  ERROR 
RETURN 

! 

:  QDFORM  - 

'  Compute  X  *  B  and  store  in 

'  that  X  is  uppe  r 

{ 

CONTINUE 

DO  4  1=1 ,N 

DO  4  J=1 .P 

triangular 

13  =  ( J-l  )*N 

Y  (13)  =  0.0 

DO  4  K=  I  .N 

+  I 

11  =  (K-l )*N 

+  I 

12  =  ( J-l  )*N 

+  K 

Y  (13)  =  Y  (13) 
CONT I NUE 

+  X  (  I  1  )  * 

Array  is  not  positive  definite!) 


*  * 

*  *  Branch  if  ma  trix  C  is  referenced. 

*  * 

IF  (CFLAG)  GO  TO  7 

*  * 

*  *  t 

*  *  Compute  Y  *  Y  and  store  in  I)  taking  advantage  of  the 

*  *  s vmme  try  of  D . 


*  * 


DO  6 

1  =1  .  P 

DO  6 

J  =  1  .  P 

I  1 

=  (  I  -  1  )  *  \ 

I  2 

=  (J-l  )  *  N 

I  3 

=  (J-l  )  *  P 

a 

s 


.V 

.V 


A 


V) 

'J 


ft 


i 


Si 


£ 


IF 

k 


is 


a 


i 
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D  (13)  -  0.0 
DO  5  K-l .N 

D  (13)  -  D  (13)  +  Y  ( I 1+K)  *  Y  (I2+K) 

CONT I  NL'E 

14  -  ( I-1)*P  +  J 

D  (14)  -  D  (13) 

CONTINUE 

RETURN 

t  * 

**  Branch  depending  on  the  value  of  ADDFLG 

*  * 

CONTINUE 

IF  (ADDFLG)  GO  TO  10 


Compute  C  -  Y  *  Y  and  store  in  D  taking  advantage  of 

*  the  symmetry  of  D. 

* 

DO  9  1-1 .P 

DO  9  J-I .P 

11  -  (  I  -1  )  *N 

12  -  ( J-l  )*N 

13  -  ( J-l  )*P  +  I 

D  (13)  -  C  (13) 

DO  8  K-l .N 

D  (13)  -  D  (13)  -  Y  ( I 1+K)  *  Y  (I2+K) 

CONTINUE 

14  -  (I-1)*P  +  J 

D  (14)  -  D  (13) 

CONTINUE 

RETURN 


Compute  C  +  Y  *  Y  and  store  in  D  taking  advantage  of 
the  s  ymme  try  of  D. 


*  * * 

10 

CONTINUE 

DO  12  1-1  .P 

*  • 

:* 

DO  1 2  J -  I  .  P 

11  -  (  I  -1  )*N 

12  =  ( J-l  )  *N 

r*. 

13  -  ( J - 1 )*P  +  I 

£ 

D  (13)  -  C  (13) 

DO  1 1  K-l  N 

0(13)  -  D  ( I 3 )  +  Y 

V 

.4 

1 1 

CONT I NUE 

V 

14  -  ( 1 -1 )*P  +  J 

D  (14)  -  D  (13) 

i 

A 

A 

1  2 

CONT  I NUF. 

RETURN 

END 

Y  ( I 1+K)  »  Y  ( I2+K) 


*  ■'*  -V  r  »  * 


nr, 


V..SVW- 


.V  A 
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APPENDIX  B 

CONTRACTION  MAPPING  RESULTS 


Many  smaller  results  were  used  in  proving  the  existence  of  a  region  where  coupled  discrete¬ 
time  Riccati  equations  constitute  a  contraction  mapping.  This  appendix  is  a  collection  of  all  the 
minor  results  needed  for  the  derivation.  In  the  beginning  of  this  appendix  are  five  lemmas.  Each 
subsequent  lemma  builds  upon  the  results  of  previous  lemmas.  Next.  Theorem  3.2,  which  was 
stated  in  Subsection  3.4.2.  is  proved.  Then.  Facts  3.2  and  3.3  (also  found  in  Subsection  3.4.2)  are 
proved.  The  end  of  the  appendix  contains  a  summary  page  which  serves  as  a  quick  reference  to  key 
results.  For  ease  of  notation  the  2-subscript  on  the  Euclidean  norm.  I  I2  .  will  be  omitted  in  all  the 
lemmas. 


Lemma  1  ;  Given  R,  =  -Li.  €,>0.  i=1.2  then  NI^X))"1  -(r.Q'))-1  I  ^  (e*)2  Ifl.l2  IX-KI 

*1 

where  X  and  >'  are  any  two  positive-semidefinite  symmetric  matrices  and  r,(X  )  ^  R ,  +  B\ XB ,  . 


Proof  :  It  is  straightforward  to  calculate  : 

(r.(x))-1 -(r.cnr1  =  (/?,+  flfx*,)-1  -(/?,  +  btjbxx 
=  (I  +  R-xB\XBx)~XRrX  -  RrX  (I  +  B{YB,RrX  )_1 
=  (I  +  flf^fXfl,)-1  I*,-1  -  (I  +  RrXB\XBx)RrX  (I  +  BfYB.Rr1  H 
=  (I  +  R  B^XB XX  \  RrX  (I  +  B’YBtRrX  )  -  (I  +  R  ~XB(XB t)R  ~x 


(I  +  B(YB,R,~l  ) 


-1  y-1 


=  (i  +  RrxBTtxB,rx  RrxBTjBtRrx  -  Rrxn’\xBxRrx 


(I  +  BrxYB  XRX  ) 


-1  ^-1 


=  (i  +  /?r1-srx^,)-1  rtxb{ 


Y  -  X 


B ,R ,~X  (I  +  B[YBxR,~X  )_1  • 


(B-l.l) 


At  this  juncture,  lei  Rx  =  —I  which  implies  that  (#,)  *  =  e,  I  .  Then,  from  (B-l.l)  it  follows 


sg 

£ 


■X 


r  ' 


v\ 


8 


\s 


5 


that 


an 
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KI^X))"1  -(W))-1  I  =  Id  +  ^fX^H  e,5f  Y  -  X  B&  (I  +  BfYB^)-1  I 

=  U,)2  KI  +  c.SfXfl,)-1  B[  |r  -  X  B.Cl  +  eiBfYB,)-1! 

<  (€i>2  1(1  +  e.jSfXS,)-1  I  I  B^  IY  -XI  I(l  +  €i5fy5i)“1l  . 

Since.  I  (I  +  e,B  [XB j)-1  I  =  -  and  both  I  and  ^[XB,  are  positive-semidefinite 


(B-1.2) 


cr(I  +  €iB[XB ,) 


100 


Lemma  2  ;  Given  R.  =  —I,  e.>0.  i=1.2  then 

€j 


1¥,(X)- V,(K)I  6,15,1 


1  +  —iB'l2  IX  +KI 
2 


ix  -ri 


(B-2.1) 


where  X  and  Y  are  any  two  positive-semidefinite  symmetric  matrices  and 
*i(X  )  ^  (T,(X  ))“^5fX  .  Furthermore,  if  X  and  Y  are  confined  to  a  closed  ball  of  radius  8j  then 
the  following  bound  is  obtained  : 


I  ¥,(X )  —  ¥,(y )  I  e,  15,1-  1  +ei5il5,l2  ]  IX  -Y\  . 

Proof  :  For  this  and  subsequent  proofs  we  will  utilize  the  matrix  identity  : 


(B-2.2) 


W  X  -  Y  Z  =  j 


(w  -y)(x  +z)  +  (w  +r)(x  -  z) 


(B-2.3) 


where  W  .  X  .  Y  .  and  Z  are  arbitrary  but  compatibly  dimensioned  matrices.  Now.  notice  that 


,(r'(X))  "  a (R,  +  B[ X5,)  ^  a(5,)  +  a(5fX5,)  ^  * 

Use  (B-2.3)  and  the  triangle  inequality  to  compute 

8  *,(X)-  *,(>  )# 


(B-2.4) 


=  _8 
2 

,  85, 


(r.(x))-1  -cr.d'))-1  5,(x+r)+  (r.(x))-1  +  (r,(nH  5,(x-m 


(r,(x )H  - (r ,(r j)-1 « ■  i x+n  +  i (r,(x )H  +  (r,(r )H  i  » x  - r  i 


An  additional  application  of  the  triangle  inequality  coupled  with  the  result  of  Lemma  1  and  (B- 
2.4)  yields  : 


¥,(X  )  -  *,(}• ) 
^  II  5 ,  II 


=  €.15, 


(e,)2  8  5 ,  «2  -  #  X  +  Y  I!  ■ «  X  -  Y  II  +  2  e,  II  X  -  Y 
€  2  •  II  X  +  F  II 


1  +  -115, 


II  X  —  >  II  . 

Hence.  (B-2.1)  is  verified.  Finally,  if  X  and  >’  are  contained  in  a  ball  of  radius  8,  .  then  II  X  II  *5  8, 
and  ii  >  II  ^  8,  Thus. 


i 

1 


I 
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p  v  ww  r 


i 


y. 


5 


NT, 

s 


v v  vxav*  -’v  v-  .  .  .-«.  -  v  v 


i  *,(x )  -  *,(y )  i  ^  €j  1 5 ,  i 


i  +  1 5 ,  i2  i  x  +  y  i 


IX  -  H 


^  €,15,1 


i  +  j.15,12  ix  i  +  iy 

2 


IX  -Y 


«$  €,15,1 


1  +  y  I5,|2.  [8,  +  Sj 


ix  -y 


=  e,15, 


1  +  €,8,15,12 


ix  - y 


'«  ^ k*’./ •  •’.  -' .  •  . 


Lemma  3  :  Given  5,  *  —I.  €j>0.  i=  1 .2  and  any  two  positive-semidefinite  symmetric  matrices  X  t 

and  X2  .  then  whenever  €,  <  - - - holds  simultaneously  for  i= 1 .2  ,  it  follows  that 

I  5, 12  ■  I  Xj  I 


I  E,(X,.X2)  I  ^ 


1  -  €^2  1 5  2 12  •  1 5  2  I2  I  X ,  II  X  2  I 


1  -  \,(X,)X2(X2) 


where 


S1(Xl.X2)  *  I-^^XOBMX^Bl 
~2(X,.  X2)  *  I  —  ^r2(X2)5j^i(X1)52  .  and 
Xi(X)  *  €j  I  5 , 12  ■  I  X  I . 


(B-3.1) 


(B-3.2a) 


(B-3.2b) 


(B-3.3) 


Furthermore,  if  X,  is  confined  to  a  closed  ball  of  radius  8,  .  i=1.2.  then  whenever 

€.  <  - - - the  following  bound  is  obtained  : 

x/2  8, 1  5,  I2 


S,(X,.X2)  I  <  2 


(B-3.4) 


Proof  :  For  this  proof,  we  use  (B-2.4)  of  Lemma  2  and  the  definition  of  ^,(X  )  to  conclude  that 


l*i(X)l  <  €,115,1  IX  I  . 


Define 


n,(X,.X2)  t  *i(Xl)B2*£X2)Bl.  and 

fl2(X,.X2)  J  ¥2(X2)5,¥j(X,)52  . 

Then.  II  |  5,(X,.  X2)  I'1  I  =  ill-  n,(X,.  X2)  l"1 


(B-3.5) 


(B-3.6a) 


(B-3.6b) 


Note  that  if  I  fi,(X,.X2)I  <  1.  then 


a(l-  fi,(X,.  X2))  ^  g:(I)  -  oTfi.CX,.  X2))  .  Suppose  that  li  fi,(X X2)  II  <  1  .  Hence. 

ii  1 5,(x„ x2)  I-1 1  =  i  1 1  —  n,(x,.  x2)  I-1  ii 


jt  /V  o  ) 


\  V.  '* '  j*.  /  V  .  ,*J*  '!•  ^  "k*\* ■  •_  V / \ * •*_*•**** .  •'  * 


1  -  I 

Next,  for  both  i-1  and  i-2  : 

ini(X,.X2)l  ^  IBil-IB2bl^1(Xl)l-l^2(X2)l 

4  IBtIIB2l-  |  €,  I  ^ ,  I  - 1  X ,  I  j  ■  |  €2  I  5  2  I  -  I  X  2  I  j 
=  €2  I  2?1  I2  •  I  #2  I2  ■  I  Xt  I  - 1  X2  I  . 

Define  1J  ^  J I  B { 1 2  •  I  X;  I  j  ,  i=l 2.  So.  whenever  fij  <  holds  simultaneously  for  all  i 
}  .  then  it  follows  that  I  flj(X j.  X2)  I  <  1  .  Consequently,  for  all  €;  6  [0.1J) .  i-1. 2 

i  Is^x j.  x2)  |_1 1  4  - ! - 

1  i-i  n,(x,.X2)i 

<  - 1 _ . 

1  -  €1€2  I  5  1  I2  Ifljl2  -IX,  I  IX2I 

In  view  of  (B-3.3)  . 

_ 1 _  =  _ 1 _ 

1  -  \1(X1)X2(^2)  1  -  €1€,I5,«2  15,  I2  IX,  I  IX,  I 

Therefore.  (B-3.1)  is  proved.  Finally,  take  €[  =  - 1 - .  Then.  A^Xj)  \^(X2)  ^  — 

J2  8JB,  I2  ‘  2 


(B-3.8) 
€  {  1.2 


(B-3.9) 

.  Thus. 


I  (E,(X ..  X2))-1  I  < 


1 


=  2  . 
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Lemma  4:  Given  Rt  =  — I,  €,>0.  i=l,2  and  any  four  positive-semidefinite  symmetric  matrices 

ATj  .  X2  .  T,  .  and  K2  •  such  that  I  X,  I  ^  S,  and  \Y  x  I  <  8,  ,  then  if  0  <€,<«;  ^ - - — —  holds 

8, 1 5, 12 

simultaneously  for  all  i  6  {  1.2}  .  it  follows  that : 


S,(X,.X2))  1  -  Si(Y1.Y2)~l 


<  V,  v2  (1  +  Vi) 
8,  (  1  -  Vi  v2 )2 


1  II 

Vi  v2  (1  +  v2 ) 
S2  (1  -  Vi  v2)2 


(B-4.1) 


iXi-Yii  +  — Li - 2  ix2-y2i 


where  vt  ^  €,  5,  I  B,  I2  <  1  . 


Proof  :  Using  notation  of  Lemma  3.  (B-3.6a-b) : 

S,(X,.  X2)  =  I  —  ftj(X i.  X2)  . 

Then  via  a  development  that  is  similar  to  Lemma  1.  (B-l.l).  we  obtain 


|s,(x,.x2)|  1  -  |h,U',. k2)|  1  =  |hj(x ,. x2)|  1  n,(x,.x2)  -  n,(r„r2)  [h,(k r,))' 


Now.  for  i=l 


n,(x,.  x2)  -  0,(1',.  y2)  =  ^(x^v^x^fl, 


j  (V,(X,)-  Vi(Y,))B2(92(X2)  +  V2(Y2))B,  + 
2  (^iCX,)^  ^i(T,))52(^2(X2)-^2(y'2))B, 


Likewise, 


1  (V2(x2)- V2(J'2))fl1(^l(X,)  +  *,(}', ))fi2  + 

0,(X,.X2)  -  fljdjJj)  =  j  {V2(X2)  +  V2(Y2))Bi(Vi(Xi)-Vi(Yi))B2 

Remark  :  ¥,(■)  is  always  a  function  of  the  first  variable  and  '?;>(•)  is  always  a  function  of  the  second 
variable. 


Making  use  of  the  triangle  inequality,  it  is  determined  that 


n,(X,.X,)  -  n!(T,.l'2)ll  ^  LlBx\-lB2 


¥,(x,)  -  *,(}',)  ii  •  ii  ’MX.)  +  ^2(i'2)  ii  + 
¥,(X  ,)  +  ¥,(}',)  ii  II  'MX.)  -  ^2(1';) « 
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Furthermore. 

l*,(Xt)  +  *,(ri)l  €j  I  Bj  I  ( I  Xj  I  +  I  Fj  I  J  <  2  €j  8, 1 5 , 1 


(B-4.2) 


Also. 

8 

R 

l¥j(Xi)-¥i(ri)l  < 

i  +  fLi^i2  ix.  +  rj  iXi-Xji 

I0 

<  e,l5,l 

1  +  €isii5ii2  |  ix.-r.i 

V— 

* 

=  €.15,1  [ 

i  +  Vi  |  iXi-rj  . 

*V 

Therefore,  putting  all  these  facts  together  yields  the  following  bound  : 

>■ 

A 

1  ni(X,.X2)  -  ni(T1.l'2)l  ^ 

82  (1  +  V|)  1  X 1 

€1e2l51l2  I52P  8j  ( i  +  v2)l  X2 

V 

V 
s' 

e, i/2(i  + 1/,)^,!2  ix,  —  r ,  1  +  I 

(B-4.3) 


€2v,  (1  +  v2)\B29-  \X2-Y2l 


(B-4.4) 


Note  :  The  derivation  of  this  upper  bound  for  I  Ot(X  \.X2)  —  fl.d'  j.  Y 2)  I  does  not  impose  any 
restrictions  on  the  magnitude  of  €i  or  €2  !  However,  in  order  to  bound 


Hxvx2)  1  -  s,(Ki.r2)  1 


we  must  now  invoke  the  definition  of 


?"  ^ - ? — _  .  Hence,  for  all  0  <  e,  <  ,  the  hypothesis  of  Lemma  3  is  satisfied  and  we 

8,  IS*  I2 

conclude  that 


||  (s,(X1.X2)|_1  II 


1  -  X1(X1)X2(X2) 

_ 1 _ 

1  -  €,e2l  5, 12  IS2I2  IXJ  IX2 

_ 1 _ 

1  -  €,8,15, 12  e,82l52«2 


II  p,Of|tx2)|  ||  ||  |sj(K,.y2)|  II  ■  ioi(x,.x2)  -  Oi(y,.y2) 

i ni(y1..y2)  -  atYj.Yji 
( 1  -  vt  v2  j2 

6,^(1  +  *l)ia1l2  VI1  «2*1(1  +  *2)I*2I* . 

(  1  -  vx  v2  )2  (  1  -  I/,  v2  )2 


V!  V2  (1  +  V|) 


>1  (1  ~  V,  V2  j2 


ix.-r.i  + 


Vi  v2  (1  +  j/2) 


S2  (1  —  v,  u2)2 


ix,-y,i 


Lemma  5  :  Given  Rk  =  — I.  €,>0. 1=1.2  and  any  two  positive-semidefinite  symmetric  matrices  Xx 

and  X2  .  with  I  X,  I  ^  8,  .  then  whenever  0  <  €|  <  Ij  ^ - - — r-  holds  for  all  i  €  {  1.2}  .  it 

“  8,1  B^ 

follows  that : 


\BxFx{Xv  X2)l<  -  yi(l  *  V2)  and 

1  ~  v\  v2 


\B2F2(Xx.X2)\^ 


£  v2  (1  +  vx) 


*  ‘  ‘  *  \-v1v2 

where  I  and  v,  ^  €*  8*  I B  -t  &  <  1  . 


(B-5.1) 


Proof  :  By  definition  : 


jHjfX,.  x2)|  1  V,(X,) 

I-52^2(X2) 

|h2(X,.X2)|_1  *2(X2) 

I-S^X,) 

IBxFx(Xx.X2)l  ^  £15,1  I  |h,(X,.  X2)J  1  II  ¥,(  X B292(X2)  I 

IB2F2(Xx.  X2)l  ^  £IB2I  I  |h2(X,.  X2)J~l  I  IV2(X2)1  -II-Bi+i(Xi)l 
where  £  ^  I  A  I . 

Since,  I  ¥,(X  )  I  <  6,15,11X1  <  c,  8, 1  5,1. 

IBtFtiXt.X^I  ^  £ex  8, 1  5, 12  I  |s,(X X2)  |_1  I  •  1 1  -  52¥2(X2)  I 

=  £  vx  \  |Hi(X,.X2)|_1  I  II-SjWI 

4  -111 —  1I-52¥2(X2)I 
\  —  vxv2 

<  — ill —  1 1  +is2i  1  ^2(x2)  1  ]  ^  (1  +12)  . 


T  ...  >  *  <.  v<.T*. 


Similarly, 


l52/r2(^i.  *2)l  <  ^€2S2I52I2I  S2(X!,  X2)  I  II  — 5i¥i(A:1)I 


=  (val  S2(Xl.X2)  I  II -£,¥,(*,)! 


«  _ f" i 

1  -  Vl  i/2 


<  •  I  1  +lg,l  l^1(X|)l  1  ^  ^  1,2(1  +t>l) 

1  -  Vl  v2  \  1  ‘  J  1  -  V.  v-> 


vV.v 


,,  ,. ...  ■■ -.-  ■/  ■.Iv  -. 


.•  /. 
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Next,  the  proof  of  Theorem  3.2  is  given.  The  significance  of  this  result  is  that  it  defines 
Lipschitz  constants  for  the  feedback  expressions.  The  proof  is  streamlined  by  utilizing  the  results 
of  the  previous  lemmas. 

Theorem  3.2  ;  Lipschitz  Constants  for  the  Feedbacks 
Given  (X  .Y).(Z.W)  t  A  .  then 

I  B ,  (  Fi(Z  .  W )  -  Fj(X .  Y )  )  b  ^  an  I X  -  Z  h  +  a12  I  Y  -  W  l? 

I  B2  (  F2U  .W)-  F2(X  .  Y  )  )  b  ^  a2llX  -Z  l2  +  a22IY  -W  l2 

£  *1*1*2  £*[*2*1*2  £  *1*2  *1*2  .  _  £*2*1*2 

where  an  =  - — - .  a12  = - — -  .  a2l  -  - — -  .  and  a22 - = - 

Si  o2  5,  o2 

Proof  :  Consider  the  first  inequality  dropping  all  2-subscripts.  By  definition  : 

85,  (  FXU.  W)-Fj(X.y)  )l 

=  II  Bx  [=,(Z.  W)|“V,(Z)[i-52^2(W)|  -  .  Y  )  |l  —  B2V2(Y  ))  A  || 

<  £IBV\2  ||  |Ht(Z.  W)|_1^,(Z  )  I-52^2(W)|  -  |s1(X.T)|~V1(X)|l-52^20'r)|  || 

The  last  term  needs  further  investigation.  Using  the  matrix  identity  (B-2.3).  we  can  write  : 

II  |h,(z.  w ) p'^^z ) (i  —  ^2^2(vv ) J  -  |~1(x.n|~Vl(x)[i-52¥2(n|  || 

[s,(Z  .  Vy)|‘V,(Z  )  -  |s!(X.  T)|_1^i(X)  •  21  —  B2(V2(Y  )  +  ¥2(W  ))  II 

2 

H,(Z.  W)  _V,(Z)  +  SI(X.1'))‘1*,(X)  ■  B2W2(Y)  -  ¥2(W))  II 

2 

.Note  that  from  (B-3.5)  and  the  fact  that  1'  .  W  €  Bj2  by  hypothesis  : 

II  I  -  B2V2(Y  )  II  ^  1  +  iT2 

and  similarly 

II-fl2V,(W)»  ^  1  +F,  . 


i 


Therefore, 


II  (h^z.  w)J"V,(z)ji-fi2^2(w)]  -  (~,(x.n|  V,(x)ji-52¥2(nJ  || 


^  (  1  +  v2) 


£,(Z,W)  ^(Z)  -  £,(X.r)|  ^(X) 


II  (h,(z . w ) )  l*,(z)  +  j^cx.n)  Vex)  ||  i*2(n  -  *2(w)i 


(i+r2)-  ||  |£x(z, w)|  V^z)  -  (£1(x.z)pVi(x) 

^  i  a  i  i 

- - -  i  y,(z )  i  +  i  Vi(x )  i  i  *2(y )  -  *2(w )  i 

2  (  1  —  vx  v2)  I 


(B-T.3.2.1) 


where  (B-4.5)  has  been  used.  Summarizing  the  results  thus  far  : 

I5j(  Fi(Z,  W)-/'1(X,y))l 

^  g(  i  +  •  ||  |e,(z.  wo|~V,(z)  -  (H1(x.r)|"Vl(x)  II  + 


iyv(z)«  +  iv,(x)i  iy2(y)  -  *2(w)i 


£,(Z  ,  W)  V,(Z)-  ~i(X.Y)  ^,(X)  II  + 


t\B,\  iB2\ 

2  (  1  -  v]  ir2) 

<  £  (  1  +  i/p  I  B ,  I 


- -_f i‘  i  *2(r )  -  *,(w)i  . 

1  -  Vl  v2 

Now.  an  additional  use  of  the  matrix  identity  (B-2.3)  on  the  first  term  of  equation  (B-T.3.2.1) 
yields  : 


|S,(Z,W))  *,(Z)  -  |s,(X.y)J  VjCX)  || 

-  II  (Ht(Z  .  w  )  |_1  -  j£,(X.y')|_1  |v,(7)  +  ¥,(X)|  + 


1-1  1-1 

=  ,(Z  .  w  )  +  =,(X  .  Y  )  ¥t(Z)  -  *,(X) 


<  6.I,  IS  Jl  II  ^E1(Z.W)^,  -  Is.CX.nl"1  II  +  - - -  -IVjCZ)  -  V,(X)I 


Ill 


Hence. 


I5j(  FX(Z.  W)-Fl(X.Y)  )l 


«*£*;(  i  +  *p-  II  [s^z.w))  1  -  [gjCx.n 


{(1  +  S)IU.I  ,  x  ,  £F,  15, 1 

- ■_ — ■  •  •  *i(z )  -  ^(x )  i  +  - — 1 — L-  i  v2(y )  -  *2(w)  i 

i  -  Fi  r2 

<  £  (  1  +  *P  - ^ I  X  -  Z  I  +  ll.Vl  IK-WI  + 

(  1  -  vY  f2Y  12 

f  (  1  +  *>P  PJ  (  1  +  Pp  |*^  *^  (  1  +  Pp 

— =-=-  -— - —  IX  -  Z  I  +  — LJ—  - -  IK  -IV  I 

1  ~  •'l  v2  1  ”  Vl  V2  ~g2 

|  (  1  +  i/p  (  1  +  Pp  PJ  v~2 

=  - - - — i-  •  - — U±=  +  l  IX-ZI  + 

(  1  ~  v[  v2  )  1  -  v\  v2 

£  PI  P2  (  1  +  i/p  u[(  \  +  Pp 

- - 3— T  -r - =-=-  +  1  ■  I  K  —  W  I 

02  (  1  —  v,  v2  )  1  ~  vtv2 

=  |  PI  (  1  +  Pp  (  1  +  Pp  ,  ^  ^  |  +  £  PI  ^  (  1  +  Pp  (  1  +  Pp  _w  t 

^1  (  1  ~  PI  Pj>  )2  "52  (  1  ~  PJ  )2 

=  iyig  ,x-zi  +  IH555.  ,y-  w, 

* i  ^2 

i  <*„  IX  -Z  I  +  «12IK  -IV  I  . 

A  completely  analogous  argument  can  be  developed  for  the  second  inequality  which  results  in  the 
appropriate  definitions  for  a21  and  a22  ■ 


E 


.\vy,:<>;yv 


.’  V/ ■  v.< 


*_•  •  .»  TJ%J»  VjbVai  IfcNfc  * 


Now.  the  proofs  of  Facts  3.2  and  3.3  are  stated. 

Fact  32.  :  Given  any  X  .Y  6  and  B,  ^  —  I .  €,>0,  then  whenever  €j  <  - - - 

€i  I  ^  i  b2 

alii  €  {1.2}.  I  *i(X  )  -  4>i(r  )  Ij  <  - 1 - IX  -Y  12  wherei/j  <  1  . 

(1  -  vt)2 


Proof  :  Under  the  given  assumptions,  it  can  be  shown  that  : 


i  <x>i(x )  -  <*>,(r )  12 


Jl  +  €,  X  BtB{ 
<  II  |l  +  e.x  B,B^ 

Since  v,  <  1  for  all  i  €  {  1.2  I  . 


-1 


X  -Y 


-1 


I  +  e,  X  B,B[ 


1-1 


|i  +  £i r|  1  ,l2 
[i+  b\  y|_1 


i 


I  X  -  Y  I, 


cr(  I  +  €,  X  B,B[) 

1 

1  -  a(e,  X  B,Bf  ) 
1 

1  -aie,  8,  B,B[) 
1 


1  -  v, 


where  <x(  •  )  =  I  •  l2  is  the  largest  singular  value  of  (  )  .  Similarly. 


I  +  €,  B,B{  Y 


-1 


1 


l-*', 

Hence,  in  view  of  (B-F.3.2.2)-(B-F.3.2.3).  (B-F.3.2.1)  becomes 


holds  for 


(B-F.3.2.1) 


(B-F.3.2.2) 


(B-F.3.2.3) 


I  •iCX  )  —  ♦i0r"  '>"2  < 


1 


-  I  X  -  Y  l2  . 


Faet  3-3;  Given  any  (X ,  Y )  €  A  *  Bjt  x  Bj,  and  *  1 1 .  et>0.  then  whenever 

0  <  €j  <  1"  (!>;)  ^  — - r-  holds  for  all  i  €  {1,2}.  it  follows  that 

*il  B  M 

P,(X.y)  =  lA^X.Y)^  $  i-L+  VJ_  ±  £*[  .  Ki>  1 

1  —  v\  v2  — 

Pztt.n  =  IA2(.X.Y)  h  ^  £  1  i  £k2  .  k2>  1 

1  ~  **1  v2  — 

where  Vt  ^  €.7),  I  B{  \£  and  £  ^  I A  ^  . 


Proof  :  By  definition. 


MjfX.K)^  =  I A  -  B2F2(X  ,Y)\2  <  1 4  Ij  +  I  B2F2(X  .Y)\2 
M3(X.y,)l3  =  I A  -  ByF  ,(X  .y)l2  <  I  >4  I2  +  lfi,F,(X.y)l2  . 


Lemma  5  of  the  Appendix  states  that 


l5,fj(X.ni2  ^  and  IB2F2(X.Y)i2  $  £ 

1  ~  vxv2  2  6  1  -  v[iT2 


where  £  Z.  I  A  l2  Hence. 


IA,(X.ni2  <  £  +  £ 


v2  (1  +  Vj) 


1  -  vxv2 


=  £ 

,  x  *2(1  +  vi) 

=  £ 

(1  —  Vt  V2)  +  V2  +  V j  v2 

l  —  tTyV^ 

9 

1  -  FiV 2 

=  £ 

1  +  v2 

9 

1  -  v\  r2 

Defining  k,  *  -  *  +  .  then  I  A  ,(X  .  Y)  ^  ^  Similarly. 

1  *'2 


IA2(X.ni2  <  £ 


1  +  V. 


=  *~2£ 


1  -  v,)/. 

Since  0  <  i/,  <  l  .  then  inf  {  =  1  and  occurs  at  =  iT2  =  0  .  Thus.  k\,  iF,  ^  1 
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Some  Important  Facts  Which  Summarize  the  Lemmas 
In  all  cases,  i  €  {  1,2  }  . 


1)  ||  r,(X )  ||  ^  €i  for  all  X  and  all  €,  >  0 


2) 


1-1 


|r,(x)J  -  |r,(n]  1  1 1  (€t)2 1  »2  ix  -yi 

for  all  X  and  Y  and  all  €,  >  0 


3)  I  ♦.(X  )  I  $  €,!/?,!  I  X  I  for  all  X  and  all  €s  >  0  . 


4)  i*i(x)  -  ¥,(nii  <  e,  115,1 


1  +  I5.82  IX  +  Y\ 


IX  -  Yi 


for  all  X  .  Y  and  all  €,  >  0  .  Moreover,  if  X  .  Y  6  B ^  then 


l^,(X)  -  *,(n«  <  €,  IB ,1  -  1+  €,8,15,12 


IX  -  YI 


II  €,llfij!  (  1  +l/i)  IX  -Kl  . 


5) 


E,(X,.X2) 


1  —  \](X  [)\i(x2)  1  —  v\  v2 


1 


whenever  €,  <  - -  where  X,  6  5* 

5,8  5.  "2 


«  •  _  » 


(B-2.4) 


(LEMMAl) 


(B-3.5) 


(LEMMA2) 


(B-4.5) 


-"k 


*1 

r. 


£ 


m 


V*. 

v\  I 


APPENDIX  C 


,v 

,  m 
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L-A-S  CONTRACTION  MAPPING  PROGRAM  RUN 

This  appendix  contains  an  L-A-S  program  to  check  for  a  contraction  mapping  and  a  sample 
run  of  the  program.  A  summary  and  discussion  of  the  data  generated  here  may  be  found  in 
subsection  3.4.3. 


I 


>  rpf  .contra 
>pro 

1  :  Linear  Quadratic  Nash  Game 

2  (rdf)«a.bl.b2.cl,c2 

3  a.bl.b2.cl.c2(out,e)= 

4  a.bl.b2.cl.c2(game)= 

5  (rdf)«=sl.rl.s2.r2 

6  rl.r2.sl.s2(out)= 

7  cl(t).sl(*).cl(*)=ql 

8  c2(t).s2(*).c2(*)=q2 

9  ql.rl.q2.r2(lq)= 

10  (rdf)*=2l.z2 

11  ql.q2(mcp)=kl.k2 

12  ; 

13  l(dsc)=one 

14  "Enter  the  total  number  of  stages  in  this  game." 

15  (inp)=ii 

16  ;  Main  Loop 

17  (stop)= 

18  a.k  1  ,k2(out.e)= 

19  kl  ,k2(lqng)=»kln.k2n 

20  k  1  n.k2n(out.e)= 

21  ;  Contraction  Mapping  Constant 

22  kl ,zl(-)*=zzl 

23  k2.z2(-)=zz2 

24  kln.kl(-)=ttl 

25  k2n.k2(-)*tt2 

26  zzl(nrm21  =xzl 

27  zz2(nrm2)=xz2 

28  tll(nrm2)=xtl 

29  tt2(nrm2)=xi2 

30  xtl  .xt2(  +  )=xt 

31  xzl  .xz2(  +  )=xz 

32  xz(inv)=xzi 

33  xt.xzi(*)=alf 

34  xzl  ,.xz2.alf(out.e)= 

35  k  1  ,k2( mcp)=zl .z2 

36  kln.k2n(mcp)=kl.k2 

37  ii.one(-)=ii 

38  u(  if  )=a 


■ 


■  W-W- 


S 

u 
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>con 

> :  Linear  Quadratic  Nash  Game 

>  (  rdf  )-*a  ,bl  ,b2  .c  1  ,c2 

Enter  name  of  the  Data  File  (DF)  for  matrix 

Opening  file  named  :  contrl.DF 

Reading  array  named  :  a 

Reading  array  named  :  bl 

Reading  array  named  :  b2 

Reading  array  named  :  cl 

Reading  array  named  :  c2 

>  a.bl  .b2,cl  ,c2(out.e)= 


4.75537e-01  4.58790e-02  -1.13295e-04 

4.58790e-02  3.44463e-01  -7.17458e-05 

-1.13295e-04  -7.17458e-05  2.50000e-01 

bl 

9.87000e+02 

1.23000e+00 

-1.01000e-03 

b2 

1.37000e+00 

-1.00000e-03 

1.00000e-05 


>  contrl 


cl 

1.00000e+00 

0.00000e+00 

0.00000e+00 

c2 

1.00000e+00 

0.00000e+00 

0.00000e+00 


0.00000e+00 

1.00000e+00 

0.00000e+00 


0.00000e+00 

1.00000e+00 

0.00000e+00 


0.00000e+00 

0.00000e+00 

l.OOOOOe+OO 


0.00000e+00 
0.00000e+00 
1 .00000e+00 


>  a  .b  1  ,b2  .c  1  .c2(  game  )= 

>(rdf)*sl.rl.s2,r2 

Enter  name  of  the  Data  File  (DF)  for  matrix  si 

Opening  file  named  :  contr2.DF 

Reading  array  named  :  si 

Reading  array  named  :  rl 

Reading  array  named  :  s2 

Reading  array  named  :  r2 

>rl.r2,sl.s2(out)= 


>  contr2 


!W! 
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rl 

11.000 


r2 

11.000 


si 

1.000  0.  0. 

0.  1.000  0. 

0.  0.  1.000 


s2 

1.000  0.  0. 

0.  1.000  0. 

0.  0.  1.000 


>cl(t)^l(*).cl(*)-ql 

#l.sl(*).cl(*)-ql 

#2,cl(*)«ql 


>c2(t).s2(*).c2(*)-q2 

#l,s2(*).c2(*)«q2 

#2.c2(*)«q2 


>ql.rl.q2.r2(lq)= 


>(  rdf  )»z  1.22 


Enter  name  of  the  Data  File  (DF)  for  matrix  2l  >contr3 
Opening  file  named  :  contr3.DF 
Reading  array  named  :  2I 
Reading  array  named  :  22 


>  q  1  .q2(  mcp)-k  1  ,k2 


>  l(dsc)«one 


>  "Enter  the  total  number  of  stages  in  this  game. 


Enter  the  total  number  of  stages  in  this  game. 


>(inp)»ii 


***  Matrix  ii  *** 

Enter  the  dimensions  of  this  matrix.  >1.1 


Enter  the  scalar  :  ii  >15 


> ;  Main  Loop 


>(stop)- 
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>con 

kl 

1.00000e+00 

O.OOOOOe+OO 

O.OOOOOe+OO 

k2 

1.00000e+00 

0.00000e+00 

0.00000e+00 


O.OOOOOe+OO 

l.OOOOOe+OO 

O.OOOOOe+OO 


O.OOOOOe+OO 

l.OOOOOe+OO 

O.OOOOOe+OO 


O.OOOOOe+OO 

O.OOOOOe+OO 

l.OOOOOe+OO 


O.OOOOOe+OO 

O.OOOOOe+OO 

l.OOOOOe+OO 


kin 

1.00205e+00  1.55971e-02  -3.14454e-05 

1.5597le-02  1.11861e+00  -4.25852e-05 

-3.14454e-05  -4.25852e-05  1.06250e+00 

k2n 

1.00205e+00  1.55969e-02  -3.14449e-05 

1.55969e-02  1.11862e+00  -4.25854e-05 

-3.14449e-05  -4.25854e-05  1.06250e+00 

xzl 

l.OOOOOe+OO 

xz2 

l.OOOOOe+OO 

alf 

1.20666e-01 


kl 

1 .00205e+00  1.55971e-02  -3.14454e-05 

1.5597le-02  1.11861e+00  -4.25852e-05 

-3.14454e-05  -4.25852e-05  1.06250e+00 


k2 

1.00205e+00 

1.55969e-02 

-3.14449e-05 


1.55969e-02  -3.14449e-05 
1.11862e+00  -4.25854e-05 
-4.25854e-05  1.06250e+00 


kin 

1.00230e+00  1.74427e-02  -3.40677e-05 

1.74427e-02  1.13265e+00  -5.02427e-05 

-3.40677e-05  -5.02427e-05  1.06641e+00 

k2n 

1 .00229e+0<)  1.74424e-02  -3.40671e-05 

1.74424e-02  1.13265e+00  -5.02429e-05 

-3.4067le-05  -5.02429e-05  1.06641e+00 


xzl 

1.20666e-01 

xz2 

1.20666e-01 

alf 

1.18329e-01 

kl 

1.00230e+00  1.74427e-02  -3.40677e-05 

1.74427e-02  1.13265e+00  -5.02427e-05 

-3.40677e-05  -5.02427e-05  1.06641e+00 

k2 

1.00229e+00  1.74424e-02  -3.40671e-05 

1.74424e-02  1.13265e+00  -5.02429e-05 

-3.40671e-05  -5.02429e-05  1.06641e+00 

kin 

1.00232e+00  1.76606e-02  -3.43086e-05 

1.76606e-02  1.13431e+00  -5.13070e-05 

-3.43086e-05  -5.13070e-05  1.06665e+00 

k2n 

1.00232e+00  1.76603e-02  -3.43080e-05 

1.76603e-02  1.13431e+00  -5.13073e-05 

-3.43080e-05  -5.13073e-05  1.06665e+00 

xzl 

1.42782e-02 

xz2 

1.42782e-02 

alf 

1.18053e-01 

kl 

1.00232e+00  1.76606e-02  -3.43086e-05 

1.76606e-02  1.13431e+00  -5.13070e-05 

-3.43086e-05  -5.13070e-05  1.06665e+00 

k2 

1.00232e+00  1.76603e-02  -3.43080e-05 

1.76603e-02  1.13431e+00  -5.13073e-05 

-3.43080e-05  -5.13073e-05  1.06665e+00 

kin 

1  (X)233e+()0  1.76863e-02  -3.43328e-05 

1.76863e-02  1.13450e+00  -5.14426e-05 

-3.43328e-05  -5.14426e-05  1.06667e+00 
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k2n 

1.00233e+00  1.76860e-02  -3.43321e-05 

1.76860e-02  1.13450e+00  -5.14429e-05 

-3.4332  le-05  -5.14429e-05  1.06667e+00 

xzl 

1.68558e-03 


xz2 

k  * 

1.68558e-03 

alf 

5 

1.18020e-01 

kl 

Cm 

1.00233e+00 

1.76863e-02 

-3.43328e-05 

r.\ 

1.76863e-02 

1.13450e+00 

-5.14426e-05 

-3.43328e-05 

-5.14426e-05 

1.06667e+00 

k2 

k 

1.002 33e+00 

1.76860e-02 

-3.43321e-05 

1.76860e-02 

1.13450e+00 

-5.14429e-05 

-3.43321e-05 

-5.14429e-05 

1.06667e+00 

kin 

H 

1.00233e+00 

1.76893e-02 

-3.43353e-05 

IP 

1.76893e-02 

1.13453e+00 

-5.14592e-05 

-3.43353e-05 

-5.14592e-05 

1.06667e+00 

*  • 

k2n 

•J 

1.00233e+00 

1.76890e-02 

-3.43347e-05 

1.76890e-02 

1.13453e+00 

-5.14595e-05 

5 

-3.43347e-05 

-5.14595e-05 

1.06667e+00 

xzl 

1.98932e-04 

xz2 

■>: 

1.98933e-04 

alf 

■r  • 

1.18016e-01 

kl 

1.00233e+00 

1.76893e-02 

-3.43353e-05 

'v 

1.76893e-02 

1.13453e+00 

-5.14592e-05 

-3.43353e-05 

-5.14592e-05 

1.06667e+00 

k2 

■ 

1.00233e+00 

1.76890e-02 

-3.43347e-05 

1.76890e-02 

1.13453e+00 

-5.14595e-05 

... 

-3.43347e-05 

-5.14595e-05 

1.06667e+00 

o 

•r«L*v-‘vf 


raw 


Vv,>  ■  {Ofv’Cv'v-' v  vo* 
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kin 

1.00233e+00  1.76897e-02  -3.43356e-0 5 

1.76897e-02  1.13453e+00  -5.14613e-05 

-3.43356e-05  -5.146l3e-05  1.06667e+00 


k2n 

1.00233e+00  1.76894e-02  -3.43350e-05 

1.76894e-02  1.13453e+00  -5.14615e-05 

-3.43350e-05  -5.14615e-05  1.06667e+00 

xzl 

2.34772 e-05 
xz2 

2.34773e-05 

alf 

1.18016e-01 


kl 

1.00233e+00  1.76897e-02 

1.76897e-02  1.13453e+00 

-3.43356e-05  -5.14613e-05 


-3.43356e-05 

-5.146l3e-05 

1.06667e+00 


k2 

1.00233e+00  1.76894e-02  -3.43350e-05 

1.76894e-02  1.13453e+00  -5.146l5e-05 

-3.43350e-05  -5.14615e-05  1.06667e+00 

kin 

1.00233e+00  1.76897e-02  -3.43357e-05 

1.76897e-02  1.13453e+00  -5.14615e-05 

-3.43357e-05  -5.14615e-05  1.06667e+00 


k2n 

1.00233e+00 

1.76894e-02 

-3.43350e-05 


1.76894e-02  -3.43350e-05 
1.13453e+00  -5.146l8e-05 
-5.146l8e-05  1.06667e+00 


xzl 

2.77068e-06 

xz2 

2.77069e-06 

alf 

1.18016e-01 

kl 

1.00233e+00  1.76897e-02  -3.43357e-05 

1.76897e-02  1.13453e+00  -5.146l5e-05 

-3.43357e-05  -5.146l5e-05  1.06667e+00 


k2 

1.00233e+00  1.76894e-02  -3.43350e-05 

1.76894e-02  1.13453e+00  -5.14618e-05 

-3.43350e-05  -5.14618e-05  1.06667e+00 

kin 

1.00233e+00  1.76897e-02  -3.43357e-05 

1.76897e-02  1.13453e+00  -5.14615e-05 

-3.43357e-05  -5.14615e-05  1.06667e+00 

k2n 

1.00233e+00  1.76895e-02  -3.43351e-05 

1.76895e-02  1.13453e+00  -5.14618e-05 

-3.4335 le-05  -5.14618e-05  1.06667e+00 

xzl 

3.26983e-07 

xz2 

3.26986e-07 

alf 

1.18016e-01 

kl 

1.00233e+00  1.76897e-02  -3.43357e-05 

1.76897e-02  *1.13453e+00  -5.146l5e-05 
-3.4335  7e-05  -5.146l5e-05  1.06667e+00 

k2 

1.00233e+00  1.76895e-02  -3.4335  le-05 

1.76895e-02  1.13453e+00  -5.14618e-05 

-3.4335  le-05  -5.14618e-05  1.06667e+00 

kin 

1.00233e+00  1.76897e-02  -3.43357e-05 

1.76897e-02  1.13453e+00  -5.14615e-05 

-3.4335  7e-05  -5.146l5e-05  1.06667e+00 

k2n 

1.00233e+00  1.76895e-02  -3.43351e-05 

1.76895e-02  1.13453e+00  -5.146l8e-05 

-3.4335  le-05  -5.14618e-05  1.06667e+00 

xzl 

3.85891e-08 

xz2 

3.85894e-08 

alf 

1.18016e-01 
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kl 

1.00233e+00 

1.76897e-02 

-3.43357e-05 

k2 

1.00233e+00 

1.76895e-02 

-3.43351e-05 

kin 

1.00233e+00 

1.76897e-02 

-3.43357e-05 

k2n 

1.00233e+00 
1.76895e-02 
-3.4335  le-05 

xzl 

4.55412e-09 

xz2 

4.554 16e-09 
alf 

1.18016e-01 

kl 

1.00233e+00 

1.76897e-02 

-3.43357e-05 

k2 

1.00233e+00 
1.76895e-02 
-3.4335  le-05 

kin 

1  00233e+00 
1.76897e-02 
-3.43357e-05 

k2n 

1  00233e+00 
1.76895e-02 
-3.4335  le-05 

xzl 

5.37457e-10 


1.76897e-02 

1.13453e+00 

-5.14615e-05 


1.76895e-02 

1.13453e+00 

-5.146l8e-05 


1.76897e-02 

1.13453e+00 

-5.14615e-05 


1.76895e-02 

1.13453e+00 

-5.14618e-05 


1.76897e-02 

1.13453e+00 

-5.14615e-05 


1.76895e-02 

1.13453e+00 

-5.14618e-05 


1.76897e-02 

1.13453e+00 

-5.14615e-05 


1.76895e-02 
1.13453e+00 
-5 . 1 46 1 8e-05 


-3.43357e-05 

-5.14615e-05 

1.06667e+00 


-3.4335  le-05 
-5.14618e-05 
1.06667e+00 


-3.4335  7e-05 
-5.14615e-05 
1.06667e+00 


-3.4335  le-05 
-5.14618e-05 
1.06667e+00 


-3.4335  7e-05 
-5.14615e-05 
1.06667e+00 


-3.4335  le-05 
-5.146l8e-05 
1.06667e+00 


-3.43357e-05 

-5.14615e-05 

1.06667e+00 


-3.4335  le-05 
-5.14618e-05 
1.06667e+00 


xz2 
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5.3746 3e-10 
alf 

1.18016e-01 

kl 

1.00233e+00 

1.76897e-02 

-3.43357e-05 

k2 

1  00233e+00 
1.76895e-02 
-3.4335  le-05 

kin 

1.00233e+00 

1.76897e-02 

-3.43357e-05 

k2n 

1.00233e+00 
1.76895e-02 
-3.4335  le-05 

xzl 

6.34284e-l  1 
xz2 

6.34291e-ll 

alf 

1.18016e-01 

kl 

1.00233e+00 

1.76897e-02 

-3.43357e-05 

k2 

1.00233e+00 
1.76895e-02 
-3.4335  le-05 

kin 

1  <X)233e+(X) 
1  76897e-02 
-3.43357e-05 

k2n 

1.00233e+00 
1  76895e-02 
-3.4335  le-05 


1.76897e-02 

1.13453e+00 

-5.146l5e-05 


1.76895e-02 

1.13453e+00 

-5.14618e-05 


1.76897e-02 

1.13453e+00 

-5.14615e-05 


1.76895e-02 

1.13453e+00 

-5.14618e-05 


1.76897e-02 

l.l3453e+00 

-5.14615e-05 


1.76895e-02 

1.13453e+00 

-5.14618e-05 


1 ,76897e-02 
1.13453e+00 
-5.14615e-05 


1.76895e-02 

l.l3453e+00 

-5.14618e-05 


-3.43357e-05 

-5.14615e-05 

1.06667e+00 


-3.4335  le-05 
-5.146l8e-05 
1.06667e+00 


-3.43357e-05 

-5.14615e-05 

1.06667e+00 


-3.4335  le-05 
-5.14618e-05 
1.06667e+00 


-3.43357e-05 

-5.14615e-05 

1.06667e+00 


-3.4335  le-05 
-5.14618e-05 
1.06667e+00 


-3.43357e-05 

-5.14615e-05 

1.06667e+(X) 


-3.4335  le-05 
-5.14618e-05 
1.06667e+00 


xzl 

7.48554e-12 


gJt  «L,.  .  Cat  I 


l,«>  *Jk  f  » 
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xz2 

7.48565e-12 

alf 

1.18014e-01 


H 

kl 

1.00233e+00 

1.76897e-02  -3.43357e-05 

1.76897e-02 

1.13453e+00  -5.146l5e-05 

>\ 

-3.4335  7e-05 

-5.14615e-05  1.06667e+00 

JL. 

k2 

'r 

> 

1.00233e+00 

1.7b895e-02  -3.4335  le-05 

v: 

1.76895e-02 

1.13453e+00  -5.14618e-05 

-3.4335  le-05 

-5.14618e-05  1.06667e+00 

*■ 

V 

kin 

v 

1.00233e+00 

1.76897e-02  -3.43357e-05 

1.76897e-02 

1.13453e+00  -5.14615e-05 

*  *, 

-3.43357e-05 

-5.14615e-05  1.0666  7e+00 

k2n 

i 

1.002 33e+00 

1.76895e-02  -3.4335  le-05 

1 

1.76895e-02 

1.13453e+00  -5.14618e-05 

-3.4335  le-05 

-5.14618e-05  1.06667e+00 

§ 

xzl 

8.83404e-13 

t 

xz2 

V 

8.83404e-13 

,v 

alf 

1.18020e-01 

kl 

1.00233e+00 

1.76897e-02  -3.4335  7e-05 

1.76897e-02 

1.13453e+00  -5.14615e-05 

-3.4335  7e-05 

-5.14615e-05  1.06667e+00 

li 

k2 

1  00233e+00 

1.76895e-02  -3.43351e-05 

v! 

1.76895e-02 

1.13453e+00  -5.14618e-05 

V, 

•  m 

-3.4335  le-05 

-5.14618e-05  1.06667e+00 

V 

kin 

1 

j  1.00233e+00 

1.76897e-02  -3.43357e-05 

1.76897e-02 

1.13453e+00  -5.14615e-05 

-  a 

• 

-3.4335  7e-05 

-5.14615e-05  1.06667e+00 

1 
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k2n 

1.00233e+00 
1.76895e-02 
-3.4335  le-05 


1.76895e-02 

1.13453c+00 

-5.14618e-05 


-3.4335  le-05 
-5.14618e-05 
1.06667e+00 


xzl 

1.04246e-13 

xz2 

1.04273e-13 

alf 

1.18080e-01 

>end 


i 
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